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Abstract 

We present a method for the solution of the Dynamic Programming Equation (DPE) 
that arises in an infinite horizon optimal control problem. The method extends to the 
discrete time version of Al'brecht's procedure for locally approximating the solution 
of the Hamilton Jacobi Bellman PDE. Assuming that the dynamics and cost are 
C r_1 (R n+m ) and C r (M. n+m ) smooth, respectively, we explicitly find expansions of the 
optimal control and cost in a term by term fashion. These finite expansions are the 
first (r — l)th and rth terms of the power series expansions of the optimal control 
and the optimal cost, respectively. 

Once the formal solutions to the DPE are found, we then prove the existence of 
smooth solutions to the DPE that has the same Taylor series expansions as the formal 
solutions. The Pontryagin Maximum Principle provides the nonlinear Hamiltonian 
dynamics with mixed direction of propagation and the conditions for a control to 
be optimal. We calculate the forward Hamiltonian dynamics, the dynamics that 
propagates forward in time. We learn the eigenstructure of the Hamiltonian matrix 
and symplectic properties which aid in finding the graph of the gradient of the optimal 
cost. Furthermore, the Local Stable Manifold Theorem, the Stokes' Theorem, and 
the Implicit Function Theorem are some of the main tools used to show the optimal 
cost and the optimal control do exist and satisfy the DPE. 

Assuming that there is a forward Hamiltonian dynamics is to assume that the 
Hamiltonian matrix is invertible. If is eigenvalue of the Hamiltonian matrix, then 
is also a closed loop eigenvalue. We consider as a possible closed loop eigenvalue 
since the optimal cost calculated term by term is true for all closed loop eigenvalues 
with magnitude less than 1. We prove that there exists a local stable manifold for 

viii 



the bidirectional nonlinear Hamiltonian dynamics. 

We also present a method for numerically solving the Hamilton Jacobi Bellman 
(HJB) PDE that arises in the infinite horizon optimal control problem. We compute 
smooth solutions of the optimal control and the optimal cost up to some degree r — 1 
and r, respectively. The first step is to approximate around using Al'brecht's local 
solutions. Then, using a Lyapunov criteria we find a point on a level set where 
we truncate the approximated solutions and begin another polynomial estimates. 
We generate the polynomial solutions in the same fashion described in the Cauchy- 
Kovalevskaya Theorem. This process is repeated over at other points as the smooth 
solutions are patched together like a circular quilt. 
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Chapter 1 



Introduction 



We solve two equations: the Hamilton- Jacobi-Bellman (HJB) PDE and the Dynamic 
Programming Equations (DPE). Both equations are associated with the infinite hori- 
zon optimal control problem of minimizing a running cost subject to a nonlinear con- 
trol dynamics. Since the optimal control problems arise ubiquitously in engineering, 
economics and biological science where models of these types surface naturally, this 
is the motivation to conjure up methods in solving them. It is actually an economist 



who first studied optimal control problems in modelling capital accumulation [28 



However, it is the advent of modern control theory that had considerable impact on 
the treatment of these problems. The classic optimal control problem, the nonlinear 
regulator problem, is what we consider in this dissertation. We seek a feedback or 
a control which minimizes the running cost under the nonlinear control system for 
the regulator problem. The optimal control maintains the dynamics close to while 
keeping the expenditure of the cost at a minimum. The Dynamic Programming tech- 
nique is used to derive the HJB PDE and the DPE from the infinite horizon optimal 
control problem. Thus, our main reason for solving the HJB PDE and DPE is to find 
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a stabilizing feedback and a Lyapunov function to confirm local asymptotic stability 
of the given nonlinear control system. The HJB PDE corresponds to the optimal con- 
trol problem with continuous dynamics and cost while the DPE is the discrete-time 
analogue. 

Around 1961 one of the earliest results was made by Albrecht U on the nonlinear 
regulator problem. He discovered the formal power series solutions to the HJB equa- 
tions while studying analytical nonlinear control systems. In 1969, Lukes proves the 
existence and uniqueness of these power series representations of the optimal control 
feedback and the optimal cost. To our knowledge there has been no extension of these 
results to the discrete-time case. So we look at the discrete-time version of the opti- 
mal control problem with the discrete-time nonlinear dynamics and cost. Using the 
method of Albrecht, we find the formal power series solutions to the DPE. Then we 
prove the existence and uniqueness of these power series solutions in Theorem [3.1.1| . 
The Pontryagin Maximum Principle (PMP) [[|, || gives the necessary condition for 
a control to be optimal. The PMP also provides the Hamiltonian dynamics satisfied 
by the optimal state and the costate trajectories. The existence of the optimal cost 
is proved in two cases. One case is when we assume that the Hamiltonian map is 
invertible. The invertibility of the Hamiltonian matrix allows us to rewrite the dy- 
namics where both the state and the costate dynamics propagate in the direction 
where time approaches infinity. The other case is when the Hamiltonian map is not 
a diffeomorphism; i.e., we have a nonlinear bidirectional Hamiltonian discrete-time 
dynamics. Finding the stable manifold for the Hamiltonian dynamics in both cases 
is key to finding the existence of the optimal cost. For invertible maps, Hartmann 
|l"3| shows the existence of a stable manifold by the method of successive approx- 
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imations on the implicit functional equation. There is the technique of using the 



Contraction Mapping Theorem on a complete space developed by Kelley [15] in 1966. 



Applications of this method can be found in the book by Carr H and the paper of 



Krener | 21fl . There is another technique by Irwin [14fl based on an application of the 



inverse function theorem on a Banach space of sequences. One of our main results is 



Theorem [4.2.1| ; we show that there exists a local stable manifold for the bidirectional 
discrete map with a hyperbolic fixed point. After a two-step process of diagonalizing 
the dynamics we apply the technique of Kelly []nj on a complete space of Lipschitz 
functions with the supremum norm. 

Computing the solutions of the HJB equations in higher dimension is major chal- 
lenge. A standard approach requires the temporal and spacial discretizations of the 
optimal control problem and then solves the corresponding nonlinear program, see 
f23fl and the appendix by Falcone in 0. Other methods for solving the HJB PDE 



are similar to those for conservation laws |26] and marching methods |29]. Although 



these standard approaches work for such equations for low dimensional systems, these 
numerical methods become infeasible for real world problems which are higher dimen- 
sional systems. Our method is a higher order approach which requires few discretiza- 
tions at each dimension. We start with Albrecht's solutions as the initial approxi- 
mations for the HJB PDE. We then truncate the computed solutions at the point 
where the optimal cost satisfies the optimality and stability constraints and begin 
new approximations at the same point. The idea is to patch together successive ap- 
proximations to obtain a larger domain for which the numerically computed optimal 



cost is still a Lyapunov function. The numerical procedure is described in |22 | 



The outline of the dissertation is as follows: In Chapter 2, we derive the DPE. 
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We first discuss the linear-quadratic regulator problem followed by the description 
of the method of finding formal power series solutions of the DPE. The coefficients 
of the polynomials are expressed in linear equations and the discrete-time algebraic 
Riccati equation (DTARE). We discuss the solvability of these equations as well. In 
Chapter 3, we prove the existence of the optimal solutions of the DPE. We show some 
properties of the Hamiltonian dynamics that are of importance in the proof of the 
existence of the optimal cost. Chapter 4 is where we prove the existence of a stable 
manifold for the noninvertible Hamiltonian map. It is actually the generalization of 
the results in Chapter 3. In Chapter 5, we summarize our numerical approach for 
solving the HJB PDE. We include numerical results from a 1-d example and describes 
the actual algorithm in more detail. 
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Chapter 2 

Formal Solution of the Dynamic 
Programming Equations in 
Discrete-Time 

We present a method for the solution of the Dynamic Programming Equations, the 
discrete-time Hamilton Jacobi Bellman PDE that arises in an infinite time optimal 
control problem. The method extends to the discrete time Al'brecht's procedure 
for locally approximating the solution. Assuming that the dynamics and cost are 
C fe_1 (R n+m ) and C k (R n+m ) smooth, we explicitly find expansions of the optimal con- 
trol and cost in a term by term fashion. These finite expansions are the first few 
terms the power series expansions of the optimal control and cost. 



CHAPTER 2. Formal Solution of the DPE 

2.1 Preliminaries 



G 



First we discuss some definitions and notions about control systems. Consider the 
n-dimensional input and output state equations 

x + = Ax + Bu 
y = Cx + Du. 

where A, B, C, and D are nx n, n x m, I x n, I x m matrices with x G M. n , y G M', 
and u G M m . We use the notation x + = Xk+i and x = Xk- 

Definition 2.1.1 A system is controllable when any inital state Xq can be driven to 
the final state xf in finite number of steps. Equivalently, if the n x nm controllability 
matrix 

C = [B, AB, . . . ,A n - 1 J B] 

has rank n, then (A, B) is a controllable pair. Otherwise (A,B) is said to be uncon- 
trollable. The pair (A,B) is stabilizable if all uncontrollable modes are asymptotically 
stable. 

If a state equation is controllable, then all eigenvalues can be assigned arbitrarily by 
introducing a state feedback. Moreover, every uncontrollable system can be diago- 
nalized into 







A c 


A 




x c 


+ 


b c 









A u 




x u 








where (A c , b c ) is controllable. If A u is stable and (A c , b c ) is controllable, then the 
system Q2.1.1D is stabilizable. As in the definition above, we refer to the eigenvalues 
of A n as the uncontrollable modes. 
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Analogously, we define observability and dectectibility, the duals of controllabilily 
and stabilizability, respectively. 

Definition 2.1.2 A sytem is observable if for any unknown xq, 3 T s.t. the input Uk, 
and output uniquely determine xq. Equivalently, if the nl x n observability matrix 

O =[C,CA,... , CA n - x Y 

has rankn. Otherwise, (A, C) is said to be unobservable. The pair (A, C) is detectable 
if the unobservable eigenvalues are stable. 

Given a control dynamics 

x + = Ax + Bu 
x(0) = xo 

where the state x G R n , the control u G M m , A, B are n x n and n x m matrices, 
respectively. We want to find u = Kx such that the system is driven to 0; i.e., the 
system 

x + = (A + BK)x 

is asymptotically stable around 0; i.e., the spectrum of A + BK lies inside the unit 
circle. We call such control a stabilizing feedback. One way to solve this stabilization 
problem is to set up an optimal control problem. This will be discuss in the next 
sections in details. 
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2.2 Discrete-Time Optimal Control Problem 

We formulate a discrete in time infinite horizon optimal control problem of minimizing 
the cost functional, 



min y~]l(x k ,Uk) 



k=0 



subject to the dynamics 



x 



+ 



f(x,u) 



x(0) = x 

where the state vector x G ffi n , the control u G M. m , and 

f(x,u) = Ax + Bu + f [2] (x,u) + f [3] (x,u) + ... 
l(x,u) = -x'Qx + x Su + -u Ru + fi 3 \x,u) + . . . 

Here we denote f^(x,u) and l^(x, u) as homogeneous polynomials in x and u of 
degree m. 

To solve the optimal control problem is to look for an optimal feedback u* = k(x) 
such that the cost functional is kept at its minimum, namely the optimal cost 7r(xo) 
of starting the system at x , while driving the dynamics to 0. 

Here are the assumptions: We assume l(x, u) is convex in x and u so that 



Q S 
S* R 



> 



and let R > 0. In addition, it is assumed that the pair (A, B) is stabilizable and the 
pair (A, Q 1 / 2 ) is detectable. 
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2.3 Derivation of the Dynamic Programming Equations 

Given that x(0) = xq the optimal value function ir(x) is defined by 



fc=0 

This value function satisfies a functional equation, called the dynamic programming 
equation. The optimal feedback k,(x) is constructed from the dynamic programming 
equation Q. First, we state the optimality principle: 

Theorem 2.3.1 Discrete-Time Optimality Principle: 



oo 





(2.3.2) 



u 



Proof: We have that 



oo 





oo 




Generalizing the optimality 



principle at the k 



-step, we have 




(2.3.3) 



it 



The optimality equation ( j2.3.2| ) is the first equation of the dynamic programming 



equations. An optimal policy u* = k(x) must satisfy 



7l(x) — 7l(f(x, U*)) — l(x, U*) = 
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if we assume convexity of the LHS of ( |2.3.3| ). We can find u* through 

d(n(x) - 7r(/(x, u)) - l(x, u) _ 
du 

which by the chain rule becomes 

-(f(x,u))-(x,u) + -(x,u) = 

Thus, 7r(x) and k(x) satisfy these equations, the Dynamic Programming Equations 
(DPE): 

ir(x) -n(f(x,u)) -l(x,u)) = (2.3.4) 
^(f(x,u))^-(x,u) + ^-(x,u) = (2.3.5) 

We now introduce a method for solving the dynamic programming equations for k(x) 
and 7r(x). 

2.4 Power Series Expansion 

Our method of solving the DPE is an extension of Al'brecht idea for continuous- 
time systems 0. We require that f(x,u) G C fe_1 (R n+m ) and l(x,u) G C k (R n+m ) are 
expressible in Taylor's form around in N £ (0). Then, all the series representation of 
the given f(x,u), l(x,u) and the unknown k(x),7t(x) are substituted into the DPE. 
With the exception of the first level, gathering the terms of the same degree in each 
equation will yield linear equations. At the first level, we obtain the discrete-time 
algebraic Riccati equation. The advantage of our process is the reduction of DPE 
which is a nonlinear system of equation into a system of one Riccati equation and 
many linear equations. We will see that the set of linear equations has a triangular 
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structure within the levels of degree, a structure that allows the system to be easily 
solvable. 

On some neighborhood around 0, N e (0) C R" +m , the dynamics and cost assume 
the following form: 

f(x,u) = Ax + Bu + f [2] (x,u) 

+f®(x,u) + ... + fW(x,u) 



l(x, u) = -x'Qx + x'Su H — u'Ru 
2 2 



+l [3] (x,u) + ... + l [d+1] (x, 



u 



where / E C^ 1 and I eC k . 

In consequence, we expect the unknowns to have power series expansions, 

tx(x) = -x'Px + n [3] (x) + ... (2.4.6) 
k{x) = Kx + k [2] (x) + ... (2.4.7) 

as well. The known matrices are A, B, /PI, ... and Q, S, R, l [3] , ... while 
P, 7r' 3 ', . . . and K, kP\ . . . are the unknowns. 

2.4.1 Special Case: Linear-Quadratic Regulator (LQR) 

The linear-quadratic regulator is an infinite-time horizon optimal control problem 
with 

f(x,u) = Ax + Bu (2.4.8) 



l(x,u) = -x'Qx + x'Su + -u'Ru; (2.4.9) 
2 2 
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i.e., the higher degree homogeneous polynomials are zero. Also, in a linear-quadratic 
regulator, we expect that the optimal cost and the optimal control will be quadratic 
and linear, respectively; i.e., 



ir(x) 
k(x) 



Kx. 



Thus, we look for P and K. After substituting all the expansions in ( |2.3.4 ) and 
collecting quadratic terms, we obtain 
1 



P - A' PA - K'B'PA - A'PBK - K'B'PBK -Q- 2SK - K'RK 



Meanwhile gathering linear terms of ( |2.3.5 ), we get 



that reduces to 



(Ax + Bu)'PB + x'S + u'R = 



K = -(B'PB + R)-\A'PB + S)'. 



x = 
(2.4.10) 

(2.4.11) 

(2.4.12) 



It follows that (|2.4.10|) can be simplified to 
1 



-x 
2 



P - A' PA + (A'PB + S){B'PB + K)-\A!PB + S)' - Q x = (2.4.13) 



by ( [2.4.12|) . Thus, ( p.4.13|) and ( p.4.11|) are the pair of equations obtained by collecting 
the quadratic terms of Q2.3.4p and the linear terms of (|2.3.5|) : 



= P-A'PA + (A'PB + S)(B'PB + R)-\A'PB + Sy-Q (2.4.14) 
K = -(B'PB + R)- l (A'PB + S)' (2.4.15) 



Equation ( 2.4.14 ) is known as the Discrete Algebraic Riccati Equation (DTARE). 
Since (B'PB + R) is positive definite, the matrix K is well-defined once P is known. 
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The next theorem gives the necessary condition for the DTARE to have a unique 
positive definite solution P. 

Theorem 2.4.1 If the pair (A, B) is stabilizable and the pair (A, Q 1 ^ 2 ) is detectable, 
then there exists a unique positive definite matrix P which satisfies the Riccati equation 
and 

\cr(A + BK)\ < 1. 
where a(A + BK) is spectrum of (A+BK). 

Moreover, the resulting feedback u = Kx is asymptotically stabilizing for the system 

x+ = {A + BK)x. 
The above theorem can be found in [0. 

2.4.2 First Level 

When the problem has nonlinear dynamics and cost with higher order terms, the first 
step is to look for the first terms of (|2.4.6|) and ( [2.4. 7|) , i.e., 



7f(x) = -x'Px + . . 

u(x) = Kx + . . . . 



The homogeneous polynomial of higher degrees of ( |2.4.6| ) and ( J2.4.71 ) are eliminated 



as the only terms gathered are the quadratic terms in ( 2.3.4 ) and the linear terms in 



Q2.3.5D when ( |2.4.(j| ) and ( |2.4.7| ) are substituted in the DPE. Then, the first level is 



exactly the linear-quadratic regulator case. Hence, we get the solutions ( [2.4. 14| ) and 

(MM- 
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2.4.3 Higher Levels: Higher Degree Terms 

In the second level, we look for the second terms of the expansions, 

7r(x) = ^x'Px + it [3] (x) + ... 
u(x) = Kx + kP\x) + . . . , 

namely ir^(x) and k^(x). Before the method is applied, we reexpressed the f(x,u) 
and l(x,u) with the first level feedback solution; i.e., 

f(x,u) = f(x,Kx + u) 
l(x,u) = l(x,Kx + u) 

These functions have power series representation through terms of 3 rd and 4 th of the 
form 

f(x,u) = (A + BK)x + Bu + f [2] (x,u) + ... 

l(x, u) = - {x'Qx + 2x'SKx + x'K'RKx) + x'Su + u'Ru + P ] (x,u) + ... 
2 



Again repeating the same process but this time grouping all the cubic terms of (2.3.4) 
gives 

n m ( x) _ n m^ A + BK)x ) = i x >( A + bk)'J® (x, o) 

+-f l2] \x, 0)P(A + BK)x + P ] (x, 0) (2.4.16) 
2 

where the terms 

x'Su + ^x'K'Ru + ^u'RKx + -x\A + BK)'PBu) + -u'B'P(A + BK)x{2A.YJ) 



CHAPTER 2. 



Formal Solution of the DPE 



15 



drops out of Q2.4.16D as it is zero by ( |2.4.12| ) and quadratic terms of ( |2.3.5| ) yields 



n [2] (x) 



[B'PB + R)- 1 ■ 



(2.4.18) 



-4— (x, 0)P(A + BK)x + B'Pp\x, 0) + B'——((A + BK)x) 

OU OJ 



dx 



x,0) 



Cancelling these terms ([2.4.17 ) is crucial because it introduces a block triangular 
structure within the levels. If appears in ( 2.4. 16|) , then the system is coupled at 
the second level. Similar condition applies for any d level. The block triangular form 
facilitates solving the system. 

The d th Level: 

We discuss the system of equations at the d th level. In the d th , the two linear 



systems are obtained by collecting the d + 1 degree terms of the [2.3.4| and d degree 
terms of the |2.3.5| . Suppose we have solved through the d — 1 th level, we incorporate 
the feedback up to the degree d — 1 into the dynamics and cost to obtain the updated 
f(x,u) and l(x,u); i.e., 

f(x,u) = f(x,Kx + K [2] (x) + K [3] (x) + ... + K [d - 1] (x) + u) 
1(X,U) = l(x,Kx + K [2] (x) + K [3] (x) + . . . + K [d ~ 1] (x) +U). 

In the d th level, the corresponding cancellation is 
x'Su + \x'K'Ru + \u'RKx + -x'(A + BK)'PBu + -u'B'P(A + BK)x =((2.4.19) 

Zi Zi Zi z 



After the substitution and collection of all the d + 1 degree terms of (|2.3.4|) and d 
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degree terms of ( |2.3.5| ), we have the equations 



n [d+1] (x) - n [d+li ((A + BK)x 



[d+l] ( 



x'(A + BK)'f [d] (x,0) 



(2.4.20) 



d-i 



+ J2[^ ] 'B' + f W (x, $)}P[B^ d - l -A + f^-^ix, 0)] 
i=2 



f [d] \x,0)P(A + BK)x 



+ n [d+1] (f) + $ d+1] (x,0) 



and 



K [d] (x) 



(B'PB + R)- 1 



du 



x,0)P(A + BK)x 



(2.4.21) 



d f)f[j] , , fl n [d+i-j] f)]\d+i] 
+ ^~(PB^ ] + + ^— (*, 0) 



i=2 



df 



du 



where ^ [d+1] (/) are terms of degree d + 1 of the homogeneous polynomials 7r' m l where 
2 < m < d + 1. These terms are generated by the input of the feedback control into 
the DPE. 



2.4.4 Solvability of the System of Equations 

We know discuss the solutions of the system of equations obtained from level > 2. 
For the case d = 2, we show that tc^ exists as a solution for (|2.4.20|) . First, the 



equation (2.4.16) is linear since 



7T [3] (x) i-> tt 13| (x) + ir [di ((A + BK)x) 



[3], 



(2.4.22) 



is linear 



Lemma 2.4.2 Given that \cr(A + BK)\ < 1. There exists a unique solution of the 
homogeneous polynomial . 
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Proof: For simplicity, we assume that (A + BK) has simple eigenvalues i.e. 

Wi(A + BK) = fjtiWi. 
where Wi G R lx ". Since the system is stable, 

\fii\ < I- (2.4.23) 
The polynomial tt^(x) has the following form: 

n n n 
«1 «2 *3 

where (wi 1 x il )(w i2 Xi 2 )(w i3 Xi 3 ) are the basis for all cubic polynomials. It follows that 



n [3] (x) -n®((A + BK)x) 



(2.4.24) 

The RHS of ( pXTSD is 

n n n 

Y Y Y d hwA w ix x ii){ w i2 x h){ w h x iz)- (2.4.25) 



ll 12 13 



Together with (|2.4.24j ) and (|2.4.25|) , we express the coefficients as 



c h,h,i3 



1 fMifH2^i3) 
where 

(1 - Hi^lHa) 7^ 

by ( p.4.23|) . Hence, there exists a unique homogeneous polynomial ir^(x). 



It follows that the linear map (|2.4.22j ) is invertible. For any level d, the homogeneous 
polynomial exists. The proof on the existence for higher order levels follows 

in the same fashion as above. 
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For any d level, the RHS of equation ( [2.4. 2ip contains a multiplicative term 



(B'PB + R)- 1 . 

The inverse matrix exists because 

(B'PB + R) 

is positive definite since R is positive definite. The rest of the terms on the RHS 
of ( j2.4.21|) are known from the previous levels. Thus, «^ exists for any level d. 



Therefore, all of the linear equations for degree 3 or higher are solvable. 
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Chapter 3 



Existence of the Local Solution of 



the DPE 



In the previous chapter, we found polynomials of degree (r — 1) and r, namely the 
optimal control and optimal cost, respectively, that satisfy the DPE. We will prove 
the existence of smooth solutions to the DPE which have the same Taylor series 
expansions as of the formal solutions previously discussed in Chapter 3. 

3.1 The Main Result 

Here we state our main result. 

Theorem 3.1.1 Suppose that the dynamics and cost are C r_1 and C r , respectively 
and the linear part of the nonlinear Hamiltonian system ( \3.4-2C\ ) is stabilizable and 
detectable around zero. Then there exists C r_1 (lR n ) optimal control and C r (lR n ) opti- 
mal cost solutions to the DPE locally around zero. 



If we let k — > oo in Theorem fl3.1.1|) , we then have the following corollary. 
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Corollary 3.1.2 Suppose that the dynamics and cost are real analytic functions on 
¥L 2n and the linear part of the nonlinear Hamiltonian system ( 3.4-2C ) is stabilizable 



and detectable around zero. Then the power series solutions converges to solution of 
the DPE locally around zero. 

3.2 PMP and Hamiltonian Dynamics 

We look at the Pontryagin Maximum Principle (PMP) as it presents the Hamilton 
difference equations satisfied by the optimal trajectories. The PMP gives a neces- 
sary condition for a control to be optimal. Along with the optimal control problem 
formulation, there is an associated nonlinear Hamiltonian, 

H(x, u, A+) = X + 'f(x, u) + l(x, u) (3.2.1) 

where A + = \k+i- The PMP states the following: 

Theorem 3.2.1 If Xk and Uk are optimal for G 0,1,2,..., then there exists A& 7^ 
for fee 0,1,2, ... such that 

x + = ^(x,u,X + ) (3.2.2) 
r) FT 

A = ^-{ x , u ,\ + ) (3.2.3) 



u 



arg min H(x, u, A + ). (3.2.4) 



Thus, the minimizer of the nonlinear Hamiltonian evaluated at the optimal x and A + 
is the optimal control u amongst all admissible controls v. Note that u*(x, A + ) G C~ x 
since / G C r ~ 1 and I G C r in ( p.2.1| ). We assume that H is convex in u to guarantee 



a unique optimal control. So the PMP gives the existence of the optimal control. 
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However, the equation ( |3.2.4| ) is an optimal control that has x and A as its independent 
variables. We must show that A is function of x to complete this part of the proof. 



3.3 Forward Hamiltonian Dynamics 



If we take the (|3.2.2| ), ( |3.2.3|) , and ( 3.2.4 ) of the PMP and linearize the equations 
around zero. The linearized sytem that we obtain is exactly the linear Hamiltonian 
system. The corresponding Hamiltonian is 

/ 1 / ili 

H(x, A + , u) = A + (Ax + Bu) + —x Qx + x Su + —u Ru 



and the system is 



x + 




X 




= e 




A 




A+ 



(3.3.5) 



where 

A-BR^S' -BR- l B' 
Q-SR- l S' A'-SR- l B' 
is the associated Hamiltonian matrix. 

Notice the opposing directions of the propagation of the state and costate dynam- 
ics in ( |3.3.5| ). This presents difficulty in studying some Hamiltonian properties, but 
it is addressed in Chapter 4. By assuming invertibility of the matrix A — SR^B', we 
can rewrite the Hamilton equations (|3.3.5| ) so that both equations propagate in the 
same direction, in particular the direction of increasing time. The forward recursion, 



(3.3.6) 



x + 




X 








A+ 




A 
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where 



(3.3.7) 



(A-BR- 1 S')-BR- 1 B'(A'-SR- 1 B')~ 1 (Q-SR- 1 S') BR' 1 B' (A' -SR' 1 B')- 1 
-(A'-SR- 1 B')- 1 (Q-SR- 1 S') (A'-SR^B')- 1 

is directly derived from the system ( |3.3.5 ). The forward recursion describes the flow 
of the state and costate as t approaches to infinity. We refer to M. F as the forward 
Hamiltonian matrix. Thus, the dynamics ( |3.4.19| ) rewritten in the forward direction 
relies on the invertibility of A' — SR~ 1 B'. 

The existence of (A' — S R -1 B') -1 has the following consequences: 

Lemma 3.3.1 If A' — SR~ 1 B' is invertible, then M F exists and is invertible. 



Proof: First, we observe in (|3~3~7l ) that (A' - SR^B')- 1 appears. Then, if (A' - 
SR^B') -1 exists so does M F . Explicit calculation gives 



(If 



F\-l 



(3.3.8) 



(A-BR- 1 S')~ 1 (A-BR- 1 S')BR- 1 B' 
-SR-^^S^iA-BR-^')- 1 (Q-SR- 1 S')(A-BR- 1 S')BR- 1 B'+(A~BR- 1 S'y 

The matrix ( |3.3.8| ) satisfies the backward hamiltonian dynamics. Thus, M B = (H F ) _1 
where 

T. T + 

(3.3.9) 

We see that H B exists since (A — BR~ 1 S') is invertible. Hence, H is invertible. a 

It follows that H F does not have zero as an eigenvalue if A' — SR'^^B' is invertible. 
Also, the Lemma (|3.3.1| ) implies that if M. F has a zero eigenvalue, then so does the 
spectrum of A — SR~ l B' . In 0, if (A, B) is stabilizable and {Q l / 2 ,A) is detectable, 
then the system is hyperbolic; i.e., none of the eigenvalues lie on the unit circle. If 



X 


= M B 


x + 






A 




A+ 
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zero is an eigenvalue, then infinity is also an eigenvalue because of the hyperbolicity 
of the system. We refer to these eigenvalues as infinite eigenvalues. This issue of 
infinite eigenvalues is discussed in Chapter 4. For the time being we will require that 
the forward Hamiltonian matrix exists and be invertible. For the rest of the chapter, 
we assume the invertibility of A' — SR^B' since the existence of H F depends on the 
inverse of A' - SR^B'. 

The forward Hamiltonian matrix has 2n eigenvalues; n are stable, i.e., these eigen- 
values lie inside the unit circle. The other n eigenvalues are unstable and are posi- 
tioned outside the unit circle. We will see in Section 3.4.3 that the eigenstructure is 
hyperbolic. 

Recall the LQR problem in Chapter 2 of minimizing 

°° 1 , 1 

min \ ^ x 'jQ x j + x'jSiij + -u'.Ruj 

3=0 

subject to the dynamics 
x + = Ax + Bu 
x = x(0). 

Lemma 3.3.2 Suppose u = Kx is the optimal control. If zero is an eigenvalue of 
(A - BR-^-S'), then zero is in a(A + BK). 

Proof: With the transformation u = Lx + v the optimal control problem is changed 



(3.3.10) 
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to the following problem: 

°° 1 1 
min V -x'jiQ + L'RL + S)xj + -v'jRvj 

j=0 

subject to the dynamics (3.3.11) 
x + = (A - BR^S^x + Bv 
x = x(0) 

where L = — i? _1 S". The transformation removes the cross term in the cost. For 
the modified problem (|3.3.11|) , we have v = Kx is the optimal control, then K = 
L + K. The transformation removes the cross term in the cost. Then, we choose the 
eigenvector xq 7^ of A — BR~ 1 S' such that its corresponding eigenvalue is zero; i.e., 

(A - BR~ 1 S')xq = (3.3.12) 

The cost starting from xq is 

tt(x )=x' (Q + L'RL + S)x . (3.3.13) 

No additional control v is exercised since the dynamics reaches zero at one time step. 
Therefore, v = Kxq = and the optimal cost is the cost (|3.3.13|) at the zero time. 
We have that 

Lx Q = Kx (3.3.14) 



for x 7^ that satisfies ( |3.3.12| ) since Kx = (L + K)x 



Returning to ( 3.3.10| ), the closed- looped spectrum is a (A + BK). Then, from 



(A + BL)x = (A + BK)x . (3.3.15) 
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Since the LHS of ( |3.3.15| ) is equal to 0, we have 

(A + BK)x = 

for xq 7^ 0. Thus, G a(A + BK). Hence, zero is a closed-loop eigenvalue. a 



By Lemma (|3.3.2|) , zero is not a closed-loop eigenvalue. This implies that zero is not 
an eigenvalue of A — BR~ 1 S' if u = Kx is assumed to be the optimal control. 

3.4 Properties of the Nonlinear Hamiltonian Dynamics 

It is imperative to look at some of the properties that the nonlinear Hamiltonian 
dynamics possesses as it is needed in the proof of the existence of the optimal cost. 
We look at its tangent dynamics, sympletic form, and eigenstructure. 

Recall that the PMP gives the nonlinear Hamilton difference equations (|3.2.2| ). 
( |3.2.3| ) satisfied by the optimal trajectories. These equations are the nonlinear Hamil- 



tonian dynamics corresponding to the optimal control problem. We calculate the 
forward nonlinear Hamiltonian dynamics from ( |3.2.2|) , (|3.2.3| ). Since the forward lin- 



ear Hamiltonian dynamics exists at zero and by the Implict Function Theorem, the 
forward nonlinear dynamics also exists in the neigborhood around 0. The forward 
nonlinear Hamiltonian dynamics, 

x + = Gx(x, A) (3.4.16) 
A+ = G 2 {x,X) (3.4.17) 

evolves in the direction of increasing time. 
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3.4.1 Tangent Dynamics 

We now introduce the idea of tangent vectors to M. = {(x,X)\x G IR n , A+ G M n )}; 
i.e., M. = M 2n . First, we denote the tangent vectors 

Sx 
SX 

The set of tangent vector to M. at (x, X) forms a vector space T Xj \M., the tangent 
space to M. at (x, A) G M. In addition, the tangent bundle of M., denoted by TA4, 
is a differential manifold where 



TM= |J T xM M. 



(3.4.18) 



The local coordinate system on TM. is 4n numbers x%, . . . , x n , X%, . . . ,X n and 

5x\, . . . , Sx n , 5X\, ... , SX n where the former is the local coordinate on M. and the 

latter is the components of the tangent vector. 

The tangent dynamics describes the flow of the tangent vector at each sequential 
point of a map. 

Definition 3.4.1 Suppose Xk+i = f{xk) with known xq. Let 8xq be a tangent vector 
at Xq. Then we define the tangent dynamics around the trajectory Xk as 

5x k+ i = ^-(x k )5x k . 

When the dynamics is linear; i.e., x k +i = Ax kl then the tangent dynamics is v k+ i = 
Avk because |£ (x k ) = A. Thus, the linear Hamiltonian dynamics ( 3.4. 19|) has 



5x+ 




5x 








5X + 




SX 



(3.4.19) 
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as it tangent dynamics. We now find the nonlinear tangent dynamics as we have 
a nonlinear system. Note the increase in difficulty in finding the tangent dynamics 
of the forward nonlinear Hamilton dynamics because of the mixed directions of the 
propagation. Linearizing the the nonlinear system at the trajectories (x, A + ) invokes 
a linear perturbation. We perturb x and A and perform the first variation on the 

dH 

o[x, A+) 

The nonlinear Hamiltonian dynamics from PMP is 

* = Q^(^ X ) 

PjTT 

A = ^(*,A+). (3.4.20) 

Suppose we replace the parameters 

x i — > x + 5x 
A + i-> A+ + 5A+. 

Then, 

x + + 5x + « (x + Sx, A + + 5X + ) 
oX + 
BM 

X + 5X « — (x + 5x,X + + 5X + ). (3.4.21) 
ox 



Subtracting (|3.4.2q ) from ( |3.4.21|) , we have 



pjTT pjTT 

6X+ ~ ^(x + 5x 1 X + + SX + )-^(x.,X + ) (3.4.22) 

pjTT pjTT 

SX « ^—( X + S X ,X + + 5X + ) - ^—(x,X + ). (3.4.23) 
ox ox 



Expanding the RHS of ( |3.4.22|) and ( g.4.23|) in Taylor series, we get 



S X = ** (l ,A +)fa ,3.4.24) 
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Equivalently, 



5x+ 
5X 



= H 5 , fc (x,A + ) 



Sx 



(3.4.25) 



(x,X + ). 



where 

H\+ x H x +\+ 
H X x H x \+ 

The subscripts of H x \+ denote partial derivatives. We look at the forward perturbed 
system. The tangent dynamics in the forward time is 



5x+ 
5\ + 



5x 
5\ 



(3.4.26) 



where 



H\+ x - H x+X +H x l x+ H xx H x+x +H xX 



~Hx\+Hxx 



To have the tangent dynamics written in forward time, we assume the invertibility of 
H x + X - When H x + X (x, A + ) is evaluated at 0, we get the matrix A' — SR^B' which we 
have assumed to be invertible. Hence, H x + X is invertible for small x and A + . 

3.4.2 The Standard Symplectic Form 

We define a nondegenerate and a bilinear symplectic two-form Vl : T( XjX )M xT( X:X )M i— > 



Q(v, w) = v'Jw and J 



/ 
-/ 



where 



Q(v , w) = —Q(w, v) 
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and 





Sx 




5x 


V = 




, w = 






5X 




T\ 



(x, A + ) G M. and (v,w) G T( X ^ X )M.. The matrix J is called the symplectic matrix. 
With the system ([3.4.26 ), it can be shown through calculations that 



K k ml k = j. 



Then, we have that 



ttT _ ttT tt-1 ttT 
n \+x n xx rL x \+ n \+\+ 



H xX+ H X+X+ 



■ H Ix H x \+ 
H xX+ 



is equal to 



where 



J 



- H\ -4— \ -4- lit 



H\+ x — H\+x+n x+x+ n xx n x + x +n xX 



H\+\+H„ 1 



H x \+Hxx 



xX+ 



(x,A + ) 



a.n ^12 
A 2 i A 22 



(3.4.27) 



A 



ii 



^12 



A. 



21 



\ 22 



H L H x \+ H x+x - H xx H xX+ H x + x +H xX+ H xx - H x+x H xX+ H xx + h xx h xX+ h x+x+ h xX+ h : 

ttT tt—T tt tt— 1 I ttT tj—1 tjT tt—T ttT tt — 1 

H xx H x\+ h * + * +H x\+ + H \+x h xX+ - h xx H xX+ H \+\+ H xX+ 

-H x \+ H \+x + H xx T + H \+\+ H xx+ H xx ~ H xX T + H x+x +H xX \H xx 

- H x\+ H *+* +H xX+ + H xX+ H X+X+ H ~ x 



xx 



-1 
A+- 



Since the matrices H x + X = H xX +, H xx , and H x + X + are symmetric, then ( |3.4.27| ) 
reduces to J. Thus, 



(3.4.28) 
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(x,A) space 



w 



w 



w 



Figure 3.1: At each sequence of the map, Q(v,w) is constant. 



Hence, the nonlinear Hamilton dynamics preserves the standard symplectic form. 
Suppose that (v, w) are two tangent vectors that satisfy ( |3.4.26|) ; i.e., v + = M.f k v and 



It! 



M.f k w. Let v, w be two sequences propagating according to( p.4.26l) . Then 



Q(v + ,w + ) = v + Jw + 



v'Jw 
Q(v, w). 



In other words, for any two sequences of tangent vector under the tangent dynamics 
( p.4.25|) , the value of the Q is unchanged at every point of the map. See fig. (|3.1| ). 
As a consequence, if we evaluate fl3.4.28p at 0, then at the linear level 



u F 'm F = j. 



(3.4.29) 



where M F ' is the matrix in ( |3.4.19| ). The calculations are as follows: For simplicity 
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we let 

a = A-BR^S' 

(3 = BR- X B' 

7 = Q-SR^S'. 



Then, 



HI 



a — (3a T 7 (3a~ 



—T -T 

-a 7 a 



e 



F' 



and 

a T — 7a -1 /? —7a" 71 
a~ l [3 aT x 

We denote a~ T as the inverse of the transpose of a. Also, the matrices (3 and 7 are 
symmetric. The equation (|3.4.29|) is 



a T — 7a T (3 —7a 



a 



/ 
-/ 



a — (3a T 7 



—T -T 

■a 7 a 



a" 1 /? 

which is simplified to 

7 — 7a _1 /3a~ T 7 — 7 + ^a~ x (3a~ T ^ I + r ya~ 1 (3a~ T — r ya~ 1 (3a~ T 

—I + a~ 1 (3a~ T ^f — a~ 1 (3a~ T ^f —a~ l (3a~ T + a~ x (3a~ T 



J. 



Therefore, 



u F 'm F = j 



holds. 
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Similarly, if we take tangent vectors satisfying the linear tangent dynamics 



v + = U F v and w + = U F w, 



then 

n(v + ,w + ) = n(v,w). 

Thus, the two-form Q(v,w) is unchanged when evaluated at the tangent vectors v 
and w of every point of the map under the linear Hamiltonian dynamics. 



3.4.3 Eigenstructure 

Now with ( |3.4.29| ), the eigenstructure of M. F is as follow: 



A 



Also, note that 



J 



= J. If fi e 


a(A) 


, then so do -, /i, and 


eigenvector of A for some eigenvalue fi. 


Sx 
5X 


= V 


5x 
5X 




5x 




-SX 




5X 




5x 





Then, 





Sx 




Sx 


/j,A'J 


Sx 




-SX 


A'JA 




= J 












SX 




SX 




SX 




Sx 



fiA' 



A 



-SX 
Sx 

-SX 

Sx 



1 



-SX 

Sx 

-SX 

Sx 
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Since \i ^ 0, we have 



/ r 



-8X 



1\ 



1 




5x J 

cr (A), the eigenvalue - is also a (A). It follows 



as the eigenpair of A. Since cr(A') 



that the complex conjugates p, and 4 G cr(A) as the characteristic polynomial p(fj) = 
p{jx) implies that p(p) = 0. Thus, the eigenvalues of A come in reciprocal pairs and 
complex conjugate pairs. B 



2n eigenvalues in the complex plane. None of the eigenvalues of H F are on the unit 
circle since the system is linearly stabilizable and detectable. 

Moreover, the eigenvectors corresponding to the stable eigenvalues sitting inside 
the unit circle spans subspace E s . Similarly, the subspace E u is spanned by the 
unstable eigenvectors whose corresponding eigenvectors lie outside the unit circle. In 
Section 3.6, we show that the subspace E s is invariant. 

3.5 Local Stable Manifold 

3.5.1 Local Stable Manifold Theorem 
Theorem 3.5.1 Given the dynamics 

x k+1 = G(x k ) 

67(0) = 

Let G : U — > W l be a C r_1 (IR 2n ) map with a hyperbolic fixed point 0. Then there is a 
local stable manifold, W s (0) G C r_1 (lR 2n ) ; that is tangent to the eigenspace Eq of the 



Since M F ' 



JM F , then the theorem above verifies the claim of the partitioning of 
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Jacobian of G atO. Define 



W s (0) = {xeU\ lim G\x) = 0} 



and 



Eq = {span of eigenvectors whose corresponding eigenvalues are such that |A| < 1}. 
W^ s (0) is a smooth manifold. 



The statement of this theorem is found in [12[]. The following theorem was proven 
by Hartman [13|], but a more modern technique can be found in ||. 

Along with our dynamics Q3.4.16 ) and a critical point 0, the local stable manifold 
theorem gives the existence of a local stable manifold W s (0). The tangency of the 
linear subspace E s to W s implies that the lowest degree term of the local stable 
manifold is a quadratic homegeneous polynomial. 



3.5.2 Construction of the Local Stable Manifold 

Through the Taylor approximation technique, we explicitly construct the local stable 
manifold term by term. These calculations coincide with the power series solutions 
found in Chapter 2. Having a system that has n finite eigenvalues outside and n 
finite eigenvalues inside the unit circle, we can find a linear transformation that block 
diagonalizes our system (|3.4.19| ). We make a linear change of coordinate such that 
the dynamics ( 3.4.19Q is block diagonalized, 



Z s 




X 




= T 




Zu 




A 
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Matrix T transforms the dynamics ( |3.4.19| ) into 



4 




A s 







z s 


+ 


fs\%si Z u ) 









A u 




Zu 




fu ( Zg ; Z u J 



(3.5.30) 



where A s is the stable matrix, A u is the unstable matrix, and 

fs(z s , Zu) fl ^ (z.sj z u) + f\ ^ {z s , Zu) + . . . 

fu(z s , Zu) f u ^ (z s , Zu) + f u ^ (z s , Zu) + . . . 

where fi^{z s , z u ) and fu\z s , z u ) are nonlinear stable and unstable homogeneous poly- 
nomials of degree d. We look for the Taylor expansion of local stable manifold 

z u = <b(z s ) G C 



-^R"); thus, we seek the form 



;^) = [21 (^) + [31 (^) + --- + ^ 1] (^)- 



(3.5.31) 



The linear term 6^(z q ) is zero as stated in the local stable manifold theorem. The 



graph z u = 4>{z s ) is said to be an invariant manifold if 

z u = 4>{ z t) whenever z u = 4>(z s ); 



(3.5.32) 



i.e., the solution of ( p.5.30|) lies in z u = <fi(z s ) for all time. By invariance (|3.5.32|) , we 
have 



Zsy Z u ) = 4>{A s Zg + f s (z s ,z u )) 
A u (p{z u ) + f u (z s , <p(z s )) = <p(A s z s + f s (z s , (f>(z s ))). 



Thus, 



</>(A a z 8 + f s (z s , 4>(z s ))) = A u (j)(z s ) + f u (z s , <p{z s )). (3.5.33) 
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Taking all the second degree terms of (|3.5.33| ), we see 

A u ^\z s ) - 4> [2] (A s z s ) = -fl 2 \z s ). (3.5.34) 

The map 

^ (z s ) h- A u cp [2] (z s ) - ^ (A s z s ) (3.5.35) 

is linear. Now, we question the invertibility of the map. The quadratic term (p^(z s ) 
of ( |3.5.34| ) as the spectrum of the map is 



where 



A u v k = £v k and WiA s = 



and 



4> [2 \z s ) = c Si>Si v k {wiZ s ) jwjZs) . (3.5.36) 



The details of the proof mirror the steps in Lemma Q2.4.2 ). Moreover, the polynomial 



0l d l(x) for all 2 < d < k also solve the following form: 

A u ^ d \z s ) - (j> [d] (A s z s ) = -fW(z a ). (3.5.37) 

Thus, (j) [d] (z s ) exists for 2 < d < n since it follows the same arguments as above. 
Hence, 4>{z s ) that is C r ~ 1 (W n ) smooth has been constructed. 

3.5.3 Lagrangian Submanifold 

In consequence, if v is in the stable subspace of M.f k , then so does M.f k v. Since 

|| j I 1 i \ . . . ' l • ■. ■->..■-- 11' 4 l W I l 1 1 1 1 ll ( } l II III ^ ( \ 1111,1 



ii \k v _> o as A; — » oo, in the limit Q(v, w) = under the dynamics (|3.4.26|) . Thus, 



Q(v, w) = if the tangent vectors are restricted in the stable subspace of Hf fe . 
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Recall the two-form Q : T (XjA) M x T {xA) M i-> R in §4.4. If v is W s C then 
so does M.f k v. Since (HTy,)^ — > as fc — > oo, in the limit io) = under the map 
( p.4.26|) . Thus, the two- form restricted to W^ s (0) is zero; i.e. 



Q(v, w) = for every v and w £ TW S . 



(3.5.38) 



In addition, W s has the maximal dimension of n. Hence, W s is a Lagrangian sub- 
manifold. 

Recall that we represent W s as the graph z u = <fi(z s ). The basis of T Zs m z \W s are 
of the form 



_d_ 

dz„. 



where G i — 1, . . . , n. 



Illuminating (|3.5.38|) , we get 





l 












1 








J 


1 










dz si 




dz sj 


d4>n 

dz si 




84>n 

dz s j 



(3.5.39) 
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where 1 on the first vector is placed on the ith row and 1 on second vector is on the 
jth row. Multiplying ( j3.5.39|) out, we have the condition 

|^-* = fori,j = l,...,n. (3.5.40) 

OZ s j OZ s i 

The equation ( |3.5.40| ) implies that <fi(z s ) is closed. Then, by the Stokes 7 Theorem 
there exists ip £ C r (W a ) such that 

<j>(z 8 ) = j^(z s ) where ^(0) = (3.5.41) 

locally on Af e (0). Thus, we have shown that the Lagrangian submanifold is the gra- 
dient of ip(z s ). 



3.6 The Optimal Cost 

From the previous section, we have found that the z u = (p(z s ) is the gradient of the 
if)(zg). In this section we find that the local stable manifold is also described by 
A = (j)(x), its original coordinates. We start by finding the linear term. 
Consider the optimal control problem of minimizing 

min - y Xx'jQkXj + 2x' i S k u k + u' k Ru k ) + x' T Px T 
u 2 

3=0 

subject to the dynamics 

Xk+i = A k x k + B k u k (3.6.42) 

x(0) = Xq 

where the state vector x £ M n , the control u £ M m , and x' t Pxt is the terminal cost. 
We denote the fundamental solution matrix, 



x k 




1 2 
X X . 


. x n 


(3.6.43) 


Afe 




A 1 A 2 . 


. X n 
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of the dynamics 







X k 


= Hp 











(3.6.44) 



X(0) 



The columns 



X' 



j = l,...,n 



are n linearly independent solutions and form a basis for the space of solutions. 
Through the transformation 



A-l — PkXk, 



the time-varying discrete Riccati equation, 

Pk+i — A' k P k+ iA k + (A' k P k+ iB k + S' k )K k + Q k , 
and the time-varying discrete control feedback, 

K k = -{B' k P k+1 B k + R k )-\A' k P k+1 B k + S k )', 
are derived as shown in the proof below. 

Theorem 3.6.1 Suppose (X k ,A k ) is solution to ( 3. 6.44 ) an d X k is invertible, then 



P k = A' k X k 1 satisfies the time-varying Riccati equation. 



Proof: From the Pontryagin Maximum Principle, we have 

U k = -R k l (B' k A' k+l + S' k X k ). 
By letting A' k = P k X k and substituting the dynamics (|3.6.42|) , 



(3.6.45) 



U k = —R k B' k P k+ i(A k X k + B k U k ) — R k S' k X-, 



(3.6.46) 
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which reduces to 

(J + R k l B' k P k+1 B)U k = -R^B' h P k+1 A k X k - R^S'.X,. (3.6.47) 



Multiplying ( |3.6.47| ) with R k becomes 

U k = K k X k . (3.6.48) 

where 

K k = ~(B' k P k+1 B k + R k y\B' k P k+1 A k + S' k ), 

the time-varying control feedback. We recall the costate equation for A in the Pon- 
tryagin Maximum Principle, 

A' k = A' k A' k+1 + Q k X k + S' k U k (3.6.49) 

Substituting the dynamics ( [3.6.42D and U k = K k X k , the costate equation becomes 



A^ = A' k P k+1 A k X k + (A' k P k+l B k + S' k )K k X k + Q k X k (3.6.50) 

Since the LHS of ( |3.6.50| ) is A^ = P k X k , we obtain the time- varying Riccati equation 
as X k is cancelled out on both sides, 

P k = A' k P k+l A k + {A' k P k+1 B k + S' k )K k + Q k . (3.6.51) 

Thus, A' k = P k X k solves the time- varying Riccati equation once (X k ,A k ) is found. 



It follows that the DTARE (|2.4.14| ) can be solved. By letting k — > -oo in the 
dynamics (|3.6.44|) , each linearly independent column vectors of (X k ,A k )' converges 
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to the stable direction and thus, forming the basis for the stable subspace. If X 1 
exists, then fundamental matrix solution has the same span as 

/ 

Hence, P k — > P as k — > — oo. Therefore, we solve the DTARE, 



P = A' PA + (A'PB + S')K + Q 



as k — > — oo. Moreover, the stable linear subspace is 



E s = 



I 
P 



(3.6.52) 



Thus, A = Px is the solution of the linear Hamiltonian equations. The linear 
term solution of the nonlinear Hamiltonian system is also A = Px. We now show 
that there exists A = <j>(x) that describes the local stable manifold of the dynamics 
which has A = Px is its linear part. 

Define the mapping F : Af e (0, 0) C K 2n — ► K" by 

F(x, A) = A - Px + g{x, A) 

where g(x, A) contains only the higher order terms. Then, for every {x, A) G A/" e (0, 0) 
F(x,X) = 0. Also the differential operator §f(0) = / is invertible since §f(0) = 
because g only contains the nonlinearities. We then invoke the IFT. It follows that 
there exists A = (j)(x) where (j)(x) G U where U is some subset of 1R™ such that 
F(x,<f>(x)) = for (x,4>(x)) G A r e (0, 0). Thus, A = 4>{x) is the graph of the local 
stable manifold in the x, A coordinates. 
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3.7 Proof of the Main Result (Theorem |5mQ 



Proof: Recall the nonlinear Hamiltonian system ( ]3.4.2U[ ). Assuming invertibility of 
H\,x a t around gives the nonlinear Hamiltonian system ( |3.4.16| ) that propagates in 



the direction of increasing time. When the nonlinear Hamiltonian sytem is linearized 
around the critical point zero, the system resembles that of the linear Hamiltonian 
system (p.4.19 ). We know that the linear Hamiltonian system is hyperbolic from 



( |3.4.29|) and (l3A2p . We then invoke the Local Stable Manifold Theorem as we have 



the dynamics ( [3.4. 16|) with a hyperbolic point 0. Thus, there exists a local stable 



manifold W s such that the linear invariant subspace ( |3.6.52|) is tangent at 0. Since 



the linear and nonlinear Hamiltonian structure preserve symplectic form (|3.4.29|) and 



( |3.4.28| ), it follows that if Q is restricted on W s then Q is zero ( |3.5.3S| ); i.e., for u, v are 
the basis in T P WL\ then Q(u,v) = as k — > oo. Then, the manifold is a Lagrangian 
submanifold. The Lagrangian submanifold described by 

Z u = <f>(Za) 

satisfies fl3.5.40|) . By Stokes' Theorem there exists ip G C r {R n ) such that 

dip 

4>{z s ) = q^-( z s) where ^(0) = 0. 

We have also shown that A = cj>{x) describes the local stable manifold by invoking 
IFT. Moreover, the linear term of A = <p(x) is A = Px which solves DTARE (|2.4.14j ). 
It follows that 

X = m = ^(x). (3.7.53) 

Thus, the graph of the gradient of the optimal cost is the local stable manifold of the 
associated Hamiltonian dynamics (|3.4.16| ). 
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Lastly, we show that n(x) and k(x) solves the DPE (|2.3.4 ), ( |2.3.5| ). From the 
Pontryagin Maximum Principle and ( |3.7.53 ), we have 



dn 



Equivalent ly, 



u* = k{x) = arg min H (x, — (x) , v) . 

v OX 



dH dir 

— (x,-(x), K (x)) = 0. 



Then, and k{x) solve DPE2. Since 



A 



dH 

dx ' 



and (|3.7.53| ), we have 



dir dn , ,df , „. dl , _ 
— — — ( X ) \X XI ) ~\~ — [XU ) 

dx dx dx ' dx 



Integrating (|3.7.55|) w.r.t. x, we get 



tc(x) — vr(/(x, k(x))) — Z(x, k.(x)) = 



(3.7.54) 



(3.7.55) 



which DPE1 (gX| ). Thus, n and k solve the DPE ($ZX§ , (gXg ). Furthermore 
7r G C r and k G C r_1 since / G C r and / G C r_1 in the Hamiltonian. B 
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Chapter 4 



Local Stable Manifold Theorem for 



the Bidirectional Discrete-Time 



Dynamics 



In this chapter, we look at the case where the linear map ( |3.3.5| ) is not a diffeo- 



morphism. This is the case where zero is a closed loop eigenvalue and therefore 
the Hamiltonian matrix is not invertible. As in the previous case, we study the 
eigenstructure and symplectic properties of the mixed direction nonlinear Hamilto- 
nian dynamics. Ultimately, we generalize the Local Stable Manifold Theorem for a 
bidirectional discrete map with a hyperbolic fixed point. 

4.1 Discrete-Time Version of Gronwall's Inequalities 

We begin with a discrete-time version of Gronwall's inequality. These lemmas will be 
useful in the proof of the local existence of a stable manifold. 
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Lemma 4.1.1 (Finite Difference Form) Suppose the sequence of scalars {uj}^L Q sat- 
isfies the difference inequality 

u k +i < $u k + L, (4.1.1) 

then 

fc-i 

3=0 

The proof of the lemma above clearly follows from summing the equation ( |4.1.1[ ) 
fc-times. 

Lemma 4.1.2 (Summation Form) Suppose {£}^1 is a sequence that satisfies 

fc-i 



|&l <C 1 ^|^-|+C 2 

3=0 



and constants C\,Ci > 0, then 

k 
3=1 

Proof: Let s k = Yl'jZo 101- Then, the sequence {s} ( *L satisfies 

Sk+i < (1 + Ci)s k + C 2 

where C\,Ci > 0. By the discrete time form of Gronwall's inequality, 

fc-i 

\Sk\ 



< (i + c 1 ) fc | So | + c 2 ^(i + c 1 ) fc - 1 - 

i=o 
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It follows that 



l&l < (i + Ci)N + c 2 

< (l + Ci)[(l + Ci) fc | 5o | 



k-l 
3=0 

k 

3=1 



k-l 

C 2 X)(l + Ci) j 

j=0 



since I so 1 = 0. 



4.2 Local Stable Manifold for the Bidirectional Discrete-Time 
Dynamics 

In this section we prove the existence of a local stable manifold A = <f>(x) for the 
Hamiltonian dynamics, 

x+ A -BR- l B< 

A Q A' 

where x, A G M™ and zero is an eigenvalue of A. The nonlinear terms, F and G, are 
C k functions for k > 1 such that 





X 




F(x,X+) 






+ 






A+ 




G(x,X + ) 



(4.2.2) 



F(0,0)=0, G(0,0) = 



(4.2.3) 



dF 



■(0,0) =0, 



dG 



-(0,0) = 0. 



d(x, A) d(x, A) 

The proof of the existence of a local stable manifold requires the discussion on 
the stability of the nonlinear state dynamics. The next subsection deals with the 
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local asymptotic stability of the state dynamics. Consequently, we describe the diag- 
onalization of the bidirectional Hamilton system. Finally, we show the existence of 
A = <j>{x). 



p(y) 



4.2.1 Preliminaries 

First, we introduce a C°° cut-off function p(y) : W l — > [0, 1] such that 

1, if0<|y|<l 
0, if|y|>2 

and < p(y) < 1 otherwise. Then we define the functions 

F(*,A+;e) := F(xpQ, A+p(^)) 
G(x,A + ;e) := G(xp^),\ + P(^-)) 



(4.2.4) 



for x, A + G R™. Since the F(x, A + ) and G(x, A + ) agree with F(x, A + ; e) and G(x, A + ; e), 
respectively, for \x\, |A + | < e, it suffices to prove the existence of a stable manifold for 
some e > 0. 

Now, we show that there exists Ni > and N 2 > such that 

(4.2.5) 
(4.2.6) 



\F(x,X;e) 


-F(x,X;e)\ 


< N ± e 


\x — x\ + \X — X\ 


\G(x, A; e) 


— G(x, A; e) 


< N ie 


\x — x\ + | A — A| 



and 



OF 



(x, A; e) 



dF 



(x,A;e) 



d(x, A) d(x, A) 

x, A; e) — — — —\ x i ^ € ) 



d(x, A) 



d(x, A) 



< A^2 

< N 2 



x-x\ + \X-X\\ (4.2.7) 
x-5| + |A-A|]. (4.2.8) 
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Since p(y) and its partial derivatives are continuous functions with compact support 
there exists M > such that 



dp, s 

dy {v) 



d 2 P 



dy< 



(y) 



< M 

< M 



for all A G M n . We also choose M > large enough that 



dF 

— (x,X) < M\x\ 



dF 



(x,A) < M\X\ 



d 2 F 



dx i dX^ 



< M, i,j = \,2 



(4.2.9) 
(4.2.10) 
(4.2.11) 



because of the condition ( f4.2.3|) for |x|, |A| < 1. By the Mean Value Theorem, 



\F{x,X;e) -F(x,X;e)\ < \F(x, A; e) — F(x, A; e) + F(x, A; e) — F(x, A; e)| 

< \F(x, A; e) - F(x, A; e)| + \F{x, A; e) - F(x, A; e) 



< 



dF 

dx 



(6,A;e) 



\x — x + 



dF 
~dX 



|A-A| 



where £i is between x and x and £2 is between A and A. Similarly, 



dF 


[x,X; e) 


<9F 


dx 


dx 


dF 


(x, A; e) 


dF 


-dx 


~ ax 



< 
< 



<9 2 F 



<9x 2 

d 2 F 



(6, A) 



dXdx 



(£1, A) 



x — x| + 
x — x| + 



a 2 F 



dxdX 
d 2 F 



dX 2 



(x,6) |A-A| 
(x,6) |A-A|. 



Next we estimate for < e < | 



<9F 

^(6,6;e) 



*-(p(— )&>— )6) «-(—)— + p(— 

ax e e dy e e e 



< ^(7)11^1(^(7)7! + Wf)l) 

< M[M + l]e 
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and 



d 2 F 



dx 2 



-k-y(p(—Ku—K2) -rri— )— +p(— ) 
ox e e oy e e e 



+ 



<9F 



(p(^)6,-)6) 



dx e e 
< M(1 + 2M) + 8M 2 



<96 e J e cV l e J e 2 



for 161 <e. 
Let 



Ni = M 2 [M+1] 

N 2 = M(1 + 2M) + 8M 2 . 



It follows that 



dF 



d(x, A) 

d 2 F 



(6,6;e) 



-(6, 6; 6) 



< N 2 , i,3 = l,2 



dx^Xi 

fo r |6I> 161 < e - The inequalities above also hold 67 as well 
Thus, 



and 



dF , dF _ ~ 

'x, A; e) — — — e ) 



<9(x, A) d(x, A) 

x, A; ej — — — —\ x i K e 





\x — 


x| + |A-A| 






\x — 


x| + |A-A| 




< 


N 2 


\x — x\ + |A 


-A| 


< 


N 2 


\x — x\ + | A 


-A| 



(4.2.12) 
(4.2.13) 



(4.2.14) 
(4.2.15) 



d(x, A) ' ' d(x, A) 
Henceforth we suppress e and write F(x, A), G(x, A) for F(x, A; e), G(x, A; e). 
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4.2.2 Stability of the Nonlinear Dynamics 

The stable invariant manifold is described by A = Px for the linear bidirectional 
Hamiltonian dynamics where P is the solution to the discrete algebraic Riccati equa- 
tion. Then, we know from Chapter 3 that the linear term of the stable manifold for 
the nonlinear bidirectional Hamiltonian dynamics, A = <f>(x), is Px. Thus, the local 
stable manifold is of the form 

A = <j)(x) = Px + ip(x) (4.2.18) 

where ift(x) contains all the nonlinear terms. 

Suppose we substitute Q4.2.18|) into the state dynamics in (|4.2.2j ), then the non- 



linear state dynamics becomes 

(I + BR- 1 B'P)x + = Ax- BR~ 1 B'4)(x + ) + F(x, Px + + 
By the Matrix Inversion Lemma (|23), we have that 

(I + BR^B'P)- 1 = (I - B{B'PB + Ry^B'P). 
Then, it follows that 

x + = (A + BK)x + U(x,x + ) (4.2.19) 
x(0) = x 

where K = -(B'PB + R)- 1 B'P and U(x,x+) = (I + BR- 1 B'P)- 1 {F(x,Px+ + 

ip(x + )) ~ BR- 1 B'i){x + ))- 

The implicit equation above can be solved. Let T : A/" e (0) C M? n R™ such that 

T(x, x + ) = x + - (A + BK)x - f^{x, x + ) = 
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for x,x + G Af € (0). Then, for G A/" e (0) the Jacobian 



dT 

dx + 



(0) 



i - (i + BR-iB'pr (|£(o, 0(o))0(o) - Birt^o)) 



5-0 



because of the condition ( 4.2.3 ) and i/)(x + ) only contains nonlinear terms. Then, by 
the Implicit Function Theorem there exists ¥(x) such that 



x + = F(x) 



(4.2.20) 



is equivalent to the earlier state dynamics ( [4.2.19 ). Moreover, the linear term of F(x) 
is (A + BK)x, i.e.; 



F(x) = {A + BK)x + F^(x) 



because 



= -I[-(A + BK)) 
= A + BK. 

It follows that F^(x) contains only the nonlinear terms and thus, 



i^(0) = 



and 



dxi 



(0) =0 i = 1,... 



n. 



The linear part of ( [4.2.19 ) is 



(4.2.21) 



(4.2.22) 



x + = (A + BK)x. 



(4.2.23) 
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Since the eigenvalues of (A + BK) lie strictly inside the unit circle, the term (A + 
BK) k x — > as k — > oo. Thus, the system ( |4.2.23| ) is asymptotically stable. Also, 



it implies that there exists a unique positive definite P that satisfies the Lyapunov 
equation 

(A + BK)'P(A + BK) -P = -I. 

Now we show the stability of the nonlinear dynamics 

x k = (A + BK)x + F^(x) 
x(0) = 0. 



We must prove that 



lim^ = 0; 

x^O X 



i.e., given any e > and any ip(x) satisfying the conditions 

^(0) = (4.2.24) 
\ip(x) -ip(x)\ < l(e)\x-x\ (4.2.25) 

where /(e) — > as e — > 0, there exists 5 > such that 

\FA X ) 



< e whenever \x\ < 5. 



x 



We define ip(x) to be the nonlinear term of the stable manifold in ( ^.5.36[ ). The 



conditions ([4.2.24 ) and ( |4.2.25| ) will be necessary for the proof of the local stable 



manifold theorem. 
Recall that 



x+ = ¥(x) = (A + BK)x + F^(x). 
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Then, 

= F{x,x + ) = x + - (A + BK)x - U{x,x + ) 

= (A + BK)x + F#(x) - (A + BK)x - fo(x, (A + BK)x + F^x)) 
= F^x) -f^{x,{A + BK)x + F^x)). 

It follows that 

F^x) = {I + BR- l B'P)- l \F{x,P{{A + BK)x + F^{x)) + i){{A + BK)x + F^{x))) 
- BR~ 1 B'i{j((A + BK)x + F^(x))] . (4.2.26) 

Let Bi = ||(J + BBr 1 B'P)-% B 2 = \\BR~ l B% P = ||P|| and a = max; [ where 
Xi e a(A + BK) and |A;| < 1. We have the following from ( |4.2.26| ), 



\F^{x)\ < Bi^Vie[|x| + \P{A + BK)x + PF4x)+4>{(A + BK)x + F lP (x)))\] 
+ BA\il>{(A + BK)x + F^(x)))\. 
because of ( |4.2.14j ). Using the Lipschitz condition (|4.2.25|) for if)(x), 



\F4x)\ < [B 1 N 1 e + a(B 1 N 1 e\\P\\+M 1 N 1 el(e)+B 1 M 2 l(e))]\x\ 
+ [BiiV ie ||P|| + M 1 N 1 el(e) + B 1 B 2 /(e)] \F^,(x)\. 

Solving for \F^(x)\, 

B 1 N 1 e + agiiVi4Pj + M^dje) + B 2 Z(e)) 
1 ~ l-(B 1 N 1 e + B 1 N 1 el(e)+M 2 l(e)) 

Let 5 = BlNl ]~^£+*p^^^ e. For some e > and e > 0, we have that 



5 > 0. Then, 



, M| M 1 N 1 e + a(M 1 N 1 e\\P\\+M 1 N 1 d(e)+M 2 l(e)) 
1 * [X)l ~ l-(M l N 1 e + M 1 N l el(e)+M 2 l(e)) W 

BiiVie + a(BiiVie||P|| + B 1 A^ 1 e/(e) + B 2 Z(e)) 
1 - (BiiVie + BiiVie/(e) + B 2 /(e)) 

< e. 
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Thus, 

F^(x) = o(\x\). 

Now, we use the Lyapunov argument. Let v(x) = x'Px. Then, 

Av(x) = v(x + ) — v (x) 

= x + 'Px + — x'Px 

= [(A + BK)x - F^{x)]'P[{A + BK)x - F^x)} - x'Px 

= x'((A + BK)'P(A + BK) - P)x + 2x'(A + BK)'PF^(x) 

= -\x\ 2 + 2x'(A + BKyPF^x). 
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since 



and 



\F4x)\ < —\x\ 



\2x'(A + BK)'PFj,(x)\ < -\x\ 2 

3 



for some p > 0. Thus, 



Av(x) 



x 



< 0. 



Therefore, the nonlinear dynamics is locally asymptotically stable uniform for all 

tj) e x. 



4.2.3 Diagonalization of the Hamiltonian Matrix 

Recall the bidirectional nonlinear dynamics in (|4.2.2|) 



x + 




A 


-BR- X B' 




X 




F(x, 


A + ) 














+ 




A 




Q 


A' 




A+ 




G(x, 
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where x, A G M n and zero is an eigenvalue of A. The nonlinear terms, F and 67, are 
C k functions for k > 1 such that 



F(0,0)=0, G(0,0) = 



dF 



d(x,X) 



(0,0) =0, 



dG 



d(x, A) 



(0,0) = 0. 



By substituting 



(4.2.27) 



(4.2.28) 



A = Px + ip(x) 

into the state dynamics above (|4.2.2| ), we get 

x + = (A + BK)x + f i ,(x,x + ) 

where U(x, x+) = (I + BR^ 1 B' P)~ 1 (F(x, Px + + ip(x+)) - BR- 1 B'i/j(x + )). 

As we substitute Q4.2.27 ) and (|4.2.28|) into the costate dynamics in ( 4.2.2 ), we 
also add = (BK)'X + — (BK)'X + . Then, the costate dynamics becomes 

A = (A + BK)'\ + + Qx + g^(x, x + ). 

where Q = Q-K' B' P(A+BK) and g^(x,x + ) = G(x, Px + +ip(x + ))+K'B'(-tp(x + )- 
Pf^(x,x + ). 

Thus, the substitution of 

A = Px + vf)(x) 

into the dynamics (|4.2.2| ) results in a new nonlinear dynamics 

f.i.fnr Trfr^ 

(4.2.29) 

A Q (A + BK)' J [A+ J [ g i ,(x,x+) 

where 



x + 




A + BK 




X 
















+ 


A 




Q (A + BK)' 




A+ 




g f (x,x + ) 



x 



(I + BR- l B'P)- l {F{x, ij(x + )) - BR- l B'^(x + )) 
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and 

9i ,(x,x + ) = G(x^(x + )) + K'B'(-^(x + )-Pf^x,x + ). 

The nonlinear terms and g$ are C k functions for k > 1 such that 

U(Q,Q)=Q, ^(0,0) = (4.2.30) 
d f* (0,0)=0, — (0,0) = 0. (4.2.31) 



d(x,x + ) ' d(x,x 
because of (|4.2.3| ), ip(x) only contains nonlinear terms and 

Now we introduce the 2 coordinate by the transformation 

A = z + Sx (4.2.32) 

for some matrix S to block diagonalize the block lower triangular Hamiltonian matrix 
in ( |4.2.29| ). By substitution, the system ( |4.2.29| ) becomes 



x + = (A + BK)x + U(x,x + ) 
z = (A + BK)'z + + (A + BKyS(A + BK)x-Sx + Qx + h^(x,x + ) 

where 

hip(x,x + ) = (A + BK)'Sf i() (x,x + ) + g^(x,x + ). 
Observe from the z dynamics above that the terms 

(A + BK)'S{A + BK)x -Sx + Qx = 
and recall that Q = Q - K'B'P(A + BK). Indeed, 

- S + A'S(A + BK) + K'B'S(A + BK) = -Q + K'B'P(A + BK). (4.2.33) 
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We know that 

- S + A'S{A + BK) = -Q, (4.2.34) 



is the discrete-time algebraic Riccati equation (DTARE). Subtracting (|4.2.34j ) from 
( |4.2.33D , we have 



K'B'S{A + BK) = K'B'P(A + BK). 

Thus, 

S = P 

and S satisfies the DTARE. Therefore, we have a diagonalized system 



(4.2.35) 



x + 




A + BK 




X 




U(x,x + ) 












+ 


z 




(A + BK)' 








g i ,(x,x + ) 



(4.2.36) 



where 



U(x,x + ) = (I + BR~ B P)~ (F(x, tp(x )) — BR~ B ip(x )) 



and 



h^(x, x + ) = (A + BK)'Sfy(x, x + ) + g f (x, x + ). 



The nonlinear terms L, and h.i, are C k functions for k > 1 such that 



d(x, x 

because of (142T30D and (MM)- 



U(0,0)=0, ^(0,0) = 
(0,0)=0, 



d U , nf *_ n dh i> ( ,0)=0. 



d(x, x 



(4.2.37) 
(4.2.38) 



CHAPTER 4. Local Stable Manifold for the Bidirectional Discrete- Time Dynamics 58 
4.2.4 The Local Stable Manifold Theorem 

Given the original dynamics ( |4.2.2| ) we look for the local stable manifold described by 
A = <f)(x). From the Theorem fl3.5.36| ) we have already proven that the linear term of 



the stable manifold for the system ( ^4.2.2| ) is Px. where P is the solution to DTARE. 



Therefore, the local stable manifold is 

\ = Px + ip(x) (4.2.39) 

where ip(x) only contains the nonlinear term of (p(x). In the two-step process of 

diagonalization of the system ( |4.2.2| ) , we introduce the z coordinate through the 



transformation 

A = z + Sx. 

Since S = P, it must be that 

z = A — Px = Px + ip(x) — Px = ip(x). 

Then, it suffices to prove existence of the local stable manifold z = if)(x) for the 
diagonalized system ( |4.2.36|) . In order to show the existence of z — ip(x), we use 



the Contraction Mapping Principle (CMP). To invoke the CMP, we will need a map 
T : X — > X that is a contraction on a complete metric space X. 



Theorem 4.2.1 Given the dynamics in (\j.2.3(\ ) with the nonlinear terms and 



are C k functions satisfying the conditions ( 4-2.371) and Q. . 2. 38j ) and a hyperbolic fixed 
point G ~R 2n , there exists a local stable manifold z = if)(x) around the fixed point 
where ifj is a C k function. 

Proof: 
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First notice that and are cut-off functions. It follows that is also a cut-off 
function. It suffices to prove the theorem for some e > since the cut-off functions 
/^(x, x + ; e) and h^(x, x + ; e) agree with fy(x,x + ) and h^(x,x + ) for \x\, \x + \ < e. 

By (|Oll )- (|4.2.17|) , (|X2| )- (gX2g ), and (gX37D - (gX38[) , there exists iVi, N 2 > 
such that 



\fy(x,y;e) - fy(x,y;e)\ < \x - x\ + \y - y\ 

\h i> {x,y\e)-h i> {x,y\e)\ < N x e \x - x\ + \y - y\ 



and 



9U 



d(x,y) 
dh^p 



{x,y;e) 



dfi 



d(x,y) 
dh^p 



d(x,y) ' ' d(x,y) 
Moreover, from the bound (|4.2.12|) we know the following 



9U 



d(x,y) 



(6, 6; e) 



< Ne 



and 



< Ne. 



(4.2.40) 
(4.2.41) 



< ^2 


\x - x\ + \y - y\ 


(4.2.42) 


< ^2 


\x-x\ + \y- y\ 


(4.2.43) 



(4.2.44) 



(4.2.45) 



d(x,y) 

for N > and |f 2 | < e. 

Henceforth we suppress e and write f<p(x, y), h^(x, y) for f^^x, y; e), h^(x, y; e). 

Let /(e) with 1(0) = and tjj G X C C°(|x| < e) where X is space of ip : E n — ► I 
such that 



^(0) = 
(x) — ip{x)\ < l(e)\x — x\ 



(4.2.46) 
(4.2.47) 
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for x, x e B e (0) C R n . We define 



sup 

\x\ <e 



X 



(4.2.48) 



Since X C C°({|x| < e}), to show X is a complete metric space it suffices to show 
that X is closed. We take a sequence {i/) n } £ X such that ip n — > ifj in C° norm. For 
large N > n, \i/j n (x) — ip(x)\ < | for all x e B e (0). Then, 

\lp(x) - 1p(x)\ < \l/}(x) - l(> n {x)\ + \lp n {x) - l/j n {x)\ + \lp n {x) - 1p(x)\ 

< i + l( e )\x-x\ + ^. 

By letting e — > 0, we have that \tp(x) — if)(x)\ < l(e)\x — x\. Thus, if) is a Lipschitz 
function. Similarly, the condition fl4.2.46|) is easily satisfied. It follows that 



1^(0) -0| < |V>(0)-Vn(0)| + |Vn(0)-0|<e. 

Thus, ip G X. Hence X is closed. Moreover, X is a complete metric space with the 
norm defined on (|4.2.48 ). 



Solving the z dynamics in (|4.2.36 ) via the variation of constants formula, we have 



jfe-i 



(A' + K'B') k -iz k + J2( A ' + K'B!) l -%{x h x l+1 ) 

1=3 



(4.2.49) 



for j < k. Let j = and k = oo, then ( 4.2.49|) changes to 

oo 

Y^(A , + K , B')%(x l ,xi +1 ) 



1=0 



(4.2.50) 



We define a mapping T : X — ► X by 

oo 

(TtP)(x ) = Y,(A + K'B') l h^ Xl ,x l+1 ) 



(4.2.51) 



1=0 
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From this fixed point equation, we look for the solution 

(Tip)(x ) = ip(x ). 

We must show Tip G X and prove T is a contraction on X. 

Suppose xo = is the initial condition. Clearly from the equation (|4.2.20| )- (|4.2.22| ), 



Xk = for all k. Together with Xk = for all k and the condition ( |4.2.38|) /i^,(0, 0) = 0, 
we have 

Tip{0) = 0. (4.2.52) 



Hence Tip satisfies the condition (|4.2.46| ). 



We now prove the Lipschitz condition ( ^.2.47| ) for Tip. 

For and the initial conditions Xq, Xq G IR n , we denote x^ = x(k,xo,ip) to 

be the solution of the state dynamics. 

x + = (A + BK)x + f^(x,x + ) 
x(0) = x . 

Similarly, for ip G X and the initial conditions Xq G W 1 , let Xk = x(k,Xo,ip) tbe the 
solution of 

x + = (A + BK)x + f i ,(x,x + ) 

x(0) = Xq. 



Recall a = max^ \Xj\ where Xj G a(A+BK) and \Xj\ < 1. Using the estimate (|4.2.40| ), 
at one-time step 

\x k +i - x~k+i\ < a\xk - x~k\ + Ni [\x k - x k \ + \xk+i - x k +i\] ■ 
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For 1 - Nie > 0, 



and recursively, 



_ . a + Nxe. _ 

Ffc+i — ^fc+i S t> — \%k ~ x k\ 

1 — N t e 



As long as 



e < 



1 — a 

2Nt '' 



then 



a + Nie 
1 - We 



< 1. 



Thus, 



\Xk ~ X k \ < \x - X Q \. 



(4.2.53) 



Using the bounds (|4.2.41|) and (|4.2.53|) 

\Til>(x )-Til>(x )\ < 



^2(A' + K'B') l (h^(x h xi+x) - h f (xi, s J+1 )) 

1=0 

oo 

< y^q' \b^{x h xi + x) - hj,{x h x l+ i)\ 

1=0 

oo 

< y^a 1 Nxe[\xi -x t \ + \x l+1 - x l+1 \] 



1=0 

OO 



< Y] al2N l e \ x -Xq\ 
1=0 

2N±e 



< 



I — a 



\x - Xq\. 



Let 



1(e) 



2N x e 
1 — a 
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Notice that /(e) — > as e — > 0. Thus, 

|Ity(x ) -T^(x )| < /(e)|x -x | 



and so (T^) satisfies the condition ([4.2 .47 ) for e > sufficiently small. 
Hence T maps from X — ► X. 
Next we show T is a contraction on X. 
We express the solutions to the state dynamics 

x k = x(k,x ,ip) 

and 

x k = x(k,x ,ifj) 

in the implicit form, 

fc-i 

x k = {A + BK) k x + Y,(A + BK) k ~ 1 ^ U{x v x ]+l ) 

3=0 

and 

fc-i 

x k = (A + BK) k x + J2( A + BKf-^Uix^ x j+1 ), 

respectively. 

We now denote Xj = x(j, x , if)) and Xj = x(j, x , if)) be the solutions to the state 

dynamics and satisfy the implicit form equations 

fc-i 

x k = (A + BK) k x + + BK) k ; \l\.(.r r x j+1 ) 

i=o 

and 

fc-i 

x k = (A + BK) k x + ^(A + BK^-if^xj, x j+1 ), 

j=0 
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respectively. 

The estimates (|4.2.4Q| )-( [4.2.41|) with the trajectories x(j, x , ip) and x(j, x , ■0) be- 



comes 



\U(x,y) -U(x,y) | < ri(e)|3/-2/|+r 2 (e)||^-^||+r 3 (e)|x-x| (4.2.54) 
\h^{x,y) - h^{x,y)\ < r x (e)|j/ - y\ + r 2 (e)\\i/j - ip\\ + r 3 (e)\x - x\ (4.2.55) 



where 



?T (e) = n 1A l(e) + n li2 e + n li3 l(e)e 
r 2 (e) = n 2 ,ie + n 2i2 e 2 
r 3 (e) = n 3 e 

and positive constants. Observe that r-j(e) — >■ as e — >■ 0. 

At one-time step, 



\x k+1 -x k+1 \ < m 2 {e)\x k - x k \ + m 3 {e)\\4> - 

where 

a + r 3 (e) 



1 - n(e) 

and 

r 2 (e) 



1 - ri(e) 

By invoking Gronwall's inequality in finite difference form ( [4.1. 1|) and assuming that 
for some small e > 

m 2 (e) < 1, 
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then 



\x k -x k \ < m 3 (e) 



fe— 1 

.3=0 



-$\\ 



1 - m 2 (e) 



< 



m 3 (e) 



1 - m 2 (e) 

With the bounds fl4.2.55|) and (|4.2.56Q , we get 



(4.2.56) 



|T^(xo)-r^(x )| < 



Y^(A' + K'B')\h^x h x l+l ) - h$(x u x l+1 )) 

1=0 

oo 

< J~] a 1 \h^{x h xi+i) - h$(xi, x l+1 ) \ 

1=0 

oo 

< [n(e)|x J+ i - x l+1 \ + r 2 (e)||<0 - $|| + r 3 (e)\x t - xj|] 



oo 

< 

«=0 



2(r 1 (e)+r 3 (e)) 



m 3 (e) 



1 - m 2 (e) 



-V +r 2 e )U- 



It follows that 



\T^(x Q )-T$(x )\ < c{e)U 



where 



c e 



2m 3 (e)(r 1 (e) +r 3 (e)) r 2 (e) 



Notice that c(e) 



(1 — m 2 (e))(l — a) 1 — a 
as e — > 0. Thus, for sufficiently small e > 

c(e) < 1. 



Hence, T is a contraction on X for e sufficiently small. Hence there exists a unique 
ip G X such that 
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4.3 Some Properties 
4.3.1 Eigenstructure 

Recall from the previous chapter the bidirectional linear Hamiltonian dynamics 

(4.3.57) 

A J I A+ 

where 



x + 




X 




= H 




A 




A+ 



m 



A -BR- l B' 
Q A' 



Definition 4.3.1 Suppose 



^ 5x ^ 

y ^A j 



( A N 



(4.3.58) 



Then we call \x the eigenvalue of the dynamics (U -3.571) . 



In the case where we do not assume as the eigenvalue of H, we found that the forward 
Hamiltonian M F is a hyperbolic system in Theorem 4.4.2. Similarly, we would like 
to show that the linear bidirectional Hamiltonian matrix H is hyperbolic; i.e., the 
eigenvalue of H lies strictly inside and outside the unit circle. 



Theorem 4.3.2 If fi is an eigenvalue of the dynamics ( 1.3.5% satisfying the relation 



( \4-.3. 58j ), then - is also an eigenvalue o/H. 



Proof: First, we decompose EI into 



/ 
/ 



Q A' 
A -BR^B' 



EL 
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Let's call 

Q A' 
A -BR- l B' 
and notice that § is symmetric. 

From the equation ( |4.3.58| ), we have the following 



which is equivalent to 



I I ^ 






°) 


1 


Sx 


\ 1*1 J 




\ 


') 


V 


SX j 









fi 





\ 


( Sx 


( Sx 






M 












V 


'J 




1° 




J 


[sxj 


[sx, 



We denote 







0^ 




f, 





\ 


H„ = 






m 










K 


1 i 








J 



Observe that 1 is an eigenvalue of H^. It follows that 

det [I — H M ] = =$> det [I — H^] = 



where 

/ 

/ 
^ 

Again, 1 is an eigenvalue of H^; i.e 



/ 

V 



/ 



5a: ^ 



^ Sx^ 1 



From above, we have 



f A 

yl 0/ 



iSx ^ 



SX 



J 



^ Sx ^ 



\ l6X I 
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which is equal to 



Thus, 



\ -Sx i 



V 1 V 



v & / 



v 5x j 



Hence, - is an eigenvalue of EL 



Note that jj, ^ 1. The eigenvalue ji —l corresponds to the trivial eigenvector 0. Also, 



the infinite eigenvalues and oo satisfy ( f4.3.58|) due to the singularity of A. 



4.3.2 Symplectic Form 



The nonlinear tangent dynamics of the bidirectional Hamiltonian system ( |3.4.20| ) that 
we derive in Chapter 3 is 



5x + 
5\ 



% fe (z, A+) 



Sx 



(4.3.59) 



where 



H 5jfc (a:,A + ) 



H\+ x H x+X + 
H X x H x \+ 

We denote H x \+, H xx , H\+ x + as partial derivatives. Recall the nondegenerate and 
bilinear symplectic two-from Q : T( Xj x)M. x T(x,\)M. \— > M, 



f2(t>,u>) = i/Jio and J 



J 
-/ 
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where 



and 



Q(v, w 


) = -n{w, 


v) 




5x 




5x 


v = 




, w = 






5X 




6X 



(x, X + ) E M and (v, w) G T {X)X) M. Also, M = T*M where x E M. We would like to 
show that under the tangent dynamics Q4.3.59 ), the two-form Q is invariant; i.e., 

Q(v,w) = tt(v + ,w + ). (4.3.60) 

Then, 



Q(v,w) = v'Jw 

/ 5x' 5x'H xx + 5X+'H xX+ ) J 
5\ + H x \+8x + 5x'H x \+5\ 



Sx 

H xx Sx + iJ :cA +5A H 



/ 



(4.3.61) 



and 



Q(v + ,w + ) = v + Jw + 

= ( 5x'H x+x + 5\+'H' x+x+ 6\+ ) J 

= —5\ + H\+ X 8x + 5x'H\+ x 5\ 
Since fl4.3.62|) and ( f4.3.61| ) are equal, then 



' H 5+X 5x + H x +\+SX + 



I 

(4.3.62) 



Q(v,w) = Q(v + ,w + ). 

Thus, for any two tangent vectors satisfying the dynamics (|4.3.59| ), the value of Q 
does not change. 
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4.4 Lagrangian Submanifold 

We have two-form 

Sl(v, w) = -5X + 'H X+ Ji + 5x'H x+ JX + . (4.4.63) 

Using the technique in Chapter 3, we have the following tangent dynamics around 
the trajectories (xj,Xj + i) 



5x 



+i = {[A + BK) + ^{x h x^ 1 ))8x :j (4.4.64) 



q x . +1 { x i> x j+i) ox j+i 
dh^p 



5Xj = (A' + K'B')SXj+i + - — —(xj,x j+ i)5x j+ i 

OXj + i 



By the Inverse Function Theorem, we can choose e > small enough so that 

dx k+1 



I — (xk,Xk+i) is invertible for \xk\, \x k +i\ < e because 



dxk+i 
since (|4.2.3C|) . 



9ft " (0,0) = / 



It follows that the tangent state dynamics is equivalent to 

5x k+ i = (i - ( Xfc , x k+1 )) ((A + BK) + ^-(x k , x k+1 ))8x k 

\ ox k+ i J \ dx k J 



and 



5x k+l = f[(l-J^L( Xi , Xi+1 )) \ (4.4.65) 



i=0 



+ BK) + ^{x i: x l+l ))5x, 
for \x k \, \x k+1 \ < e. 

. >.ave x fc -> 0, 8(a . x+) 



As we let & — > oo, we have — ► 0, d( 9 J^ +) (0) — > 0. It follows that from ( [4.4. 65| ) 
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Since 5xk — > as k — > oo, we have that 

fi(v,to) -> 

for f , t/; restricted to the tangent dynamics. Thus, the local stable manifold is a 
Lagrangian submanifold. Similarly as in Chapter 3, we have 



j dx ^ 



for i, j = 1,... ,n. (4.4.66) 



The equation (|4.4.66|) implies that ip(x) is closed. Then, by the Stokes 7 Theorem 
there exists n G C r (M n ) such that 

Bit 

M x ) = —(x) where ^(0) = (4.4.67) 
ox 

locally on some neighborhood of 0. Since z = tp(x) = A — Px, we have that 

ax 

because P is symmetric. Thus, there exists 7r G C r such that 
where 

07T . . 07f , . 9 / 1 . _ 

Hence, the local stable manifold A is the gradient of the optimal cost for the bidirec- 
tional Hamiltonian dynamics. 
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Chapter 5 



Solution of the Hamilton Jacobi 



Bellman Equations 



Our procedure for solving HJB PDE numerically is described in | 22| , but we give a 



summary of the method and some numerical results in this chapter. 

5.1 A Method for HJB 

We solve the Hamilton Jacobi Bellman (HJB) Partial Differential Equation that arises 
in many control problems. We consider the infinite horizon optimal control problem 
of minimizing the cost 



l(x,u) dt (5.1.1) 

Jt 

subject to the dynamics 

x = f(x,u) (5.1.2) 



CHAPTER 5. Solution of the Hamilton Jacobi Bellman Equations 73 
and initial condition 

x(t) = x°. (5.1.3) 

The state vector x is an n dimensional column vector, the control u is an m dimen- 
sional column vector and the dynamics f(x,u) and Lagrangian l(x,u) are assumed 
to be sufficiently smooth. 

If the minimum exists and is a smooth function 7r(a; ) of the initial condition then 
it satisfies the HJB PDE 

min <^(x)f(x,u) + l(x,u) \ = (5.1.4) 
u yox J 

and the optimal control k(x) satisfies 

k(x) = argmin { — (x)f(x, u) + l(x, u) \ — (5.1.5) 
u ydx J 

These are expressed in terms of the Hamiltonian 

H(p,x,u) — pf(x,u) + l(x,u) (5.1.6) 

where the argument p is an n dimensional row vector. The HJB PDE becomes 

dir 

= mmH(— (x),x, u) (5.1.7) 
u ox 

dn 

k(x) = aigmmH(—(x),x,u) (5.1.8) 

u OX 

We assume that the Hamiltonian H(p,x, u) is strictly convex in u for all p, x. 
Then ( |5.1.4j , |5.1.5| ) become 



and 



dir 

— (x)f(x, k(x)) + l(x, k(x)) = (5.1.9) 
ox 



— (x)-^(x, K (x)) + ^(x, k(x)) = (5.1.10) 
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These are the equations ( 5.1.9|) ,( |5.1.10 ) that we solve for the optimal cost 7r and 
optimal control k. 



5.1.1 Al'brecht's Method Revisited 

Al'brecht M solved the HJB PDE locally around zero by expanding the problem in 
a power series, 

f(x,u) = Ax + Bu + f [2] (x,u) 

+f [3] (x,u) + ... (5.1.11) 
l(x, u) = - (x'Qx + 2x'Su + u'Ru) 

+l l3] {x,u) + l [4] {x,u) + ... (5.1.12) 
ir(x) = ^x'Px + n [3] (x) 

W 4] (x) + ... (5.1.13) 

k(x) = Kx + k [2] (x) + k [3] (x) + ... (5.1.14) 

where denotes a homogeneous polynomial of degree d. As in Chapter 3, we sub- 
situte these equations into the HJB PDE (|5.1.9| , |5.1.10|) and equate the terms of like 



degree to obtain a sequence of algebraic equations for the unknowns. 

The first level is the pair of equations obtained by collecting the quadratic terms 
of ( |5.1.9| ) and the linear terms of ( |5.1.10| ). We denote the d th level as the pair of 



equations garnered from the [d+ l] th degree terms of ( |5.1.9| ) and the d th degree terms 
of ( |5X10D . 
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Here are the first set of equations: 

= A'P + PA + Q- 

(PB + S)R-\PB + Sy (5.1.15) 

K = -R-^PB + S)' (5.1.16) 

The quadratic terms of the HJB PDE reduce to the familiar Riccati equation 
( [5.1.15 ) and the linear optimal feedback ( |5.1.16|) . 



We assume A, B is stabilizable and Q, A is detectable then the Riccati equation 
has a unique positive definite solution P and the linear feedback locally exponentially 
stabilizes the closed loop system. Moreover the optimal quadratic cost is a local 
Lyapunov function for the closed loop system. 

The d th Level 

Suppose we have solved through the d — 1 th level after repeating the process d — 1 
times. It is convenient to incorporate this solution into the dynamics and the cost. 
Let K k \x) = Kx + k^(x) + k^(x) + . . . + k^(x) and define 

f(x,u) = f(x,K d - 1] (x) +u) (5.1.17) 
l(x,u) = l(x,n d ' 1] (x) +u) (5.1.18) 

These have power series expansions through terms of degree d and d + 1 of the form 

f(x,u) = (A+ BK)x + Bu + f [2] (x,u) + ...f [d] (x, u) + . . . 
l(x,u) = i (x'Qx + Ix'SKx + x'K'RKx) 
+x' Su + u'Ru + 

(x, u) + . . . + J} d+1 ^ (x,u) + ... 
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We plug these into the HJB PDE ( |5.1.9| , [5.1.1 0| ) and find that they are satisfied 
through the d — 1 level and don't involve u. Since u = «^ (x) + . . . , the d level 
equations are 

dn [d+i] 







{x)(A + BK)x 

d 

+ —^Z ( x )f [l K%, 0) + x'PBu 



dx 

d^ d+2 ~ 



i=2 



dx 



+]} d+1] (x, 0) + x'Su + -x'K'Ru 



and 



dx 
du 



x 



4 = 2 

(x, 0) + u'R 



dx du 



Because of ( p.l.lfip , u drops out of the first equation which becomes 







dx 

"'-1 d7V [d+2 



x)(A + BK)x 



i=2 



dx 

+l ld+1] (x,0). 



After this has been solved for 7r' d+1 ](x), we can solve the second for k^i 



.T 



* ^[d+2-i] ^J[i] d J[d+l] 



<9x 5m 



9m 



(5.1.19) 



(5.1.20) 



(5.1.21) 



(5.1.22) 



These equations admit a unique solution up to the smoothness of / and I since A+BK 
has all its eigenvalues in the left half plane. If / and / are real analytic, the power 



series converges to the solution of the HJB PDE locally around x = ||25|| , pO|| . The 
higher degree equations are linear; hence they are easily solvable. 
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Figure 5.1: Around each point, we generate a polynomial solution on each patch. 

Albrecht's approach develops the Taylor series expansion along a point. The 
Taylor series only converges on a small neigbhorhood and may quickly diverges outside 
that neighborhood. Also increasing the degree of the polynomial solutions does not 
guarantee a larger region where the approximated solutions are true. In our approach, 
we use Albrecht's method to generate the initial approximations around 0. Then 
the initial approximations are then improve by piecing successive approximations of 



smooth solutions around the points on the level sets, see Fig.( |5.lD . Instead we truncate 
the computed solutions at the point where the optimal cost satisfies the optimality 
and stability constraints and begin new approximations at the same point. 



5.1.2 The Lyapunov Criterion 

Let 7r d+1 l and denote the approximate cost and feedback to degrees d + 1 and d 
respectively. We assume l(x,u) > 0. We find the largest sublevel set 



7T 



x) < C 



(5.1.23) 
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such that 

dn d+i] 



(x)f(x,K d] (x)) < -(l-ei)/(x,K ld| (a;)) (5.1.24) 



dx 



-(l-e 2 )l(x,K^(x)) < -^—( x )f( x ,K d \x)) (5.1.25) 

The parameter €\ controls the rate of exponential stability of the closed loop system 
while €2 dictates the rate of the optimality of the feedback. We determine a sublevel 
set of the approximate cost on which it is an acceptable Lyapunov function and the 
approximate feedback is stabilizing and satisfying the optimality conditions. We em- 
phasize the stabilizing property of the control law rather than its optimality because 
usually optimality is only a tool to find a stabilizing feedback. Typically the goal 
is to find a control law to stabilize the system and the optimal control problem is 
formulated as a way of finding one. 

5.1.3 Power Series Expansions 

We specialize to problems with dynamics that is affme in u and with a Lagrangian 
that is quadratic in u, 

f(x,u) = g (x)+g 1 (x)u (5.1.26) 
l(x,u) = lo(x) + li(x)u + m'/ 2 (x)m (5.1.27) 

where h{x) is an invertible m x m matrix for all x. 

Assume that we have solved the HJB PDE locally around and have choosen a 



sublevel set of value c subject to the condition (|5.1.23| , |5.1.24j ). We generate a power 
series solution around a point x on the level set ir d+1 ^ (x) = c. 

We introduce some notation. Let a = (ai, . . . , a n ) be a multi- index of nonnega- 
tive integers and \a\ = Ylii a i- Let (3 = (fa, . . . , /3 n ). We say (3 < a if $ < a iy i = 
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1, . . . , n and j3 < a if (3 < a and for at least one i, $ < a^. Let = (0, 
Define the differential operator 



..,o). 



the multifactorial 



the monomial 



and the coefficient 



D a = 



d 
dxi 



"OXr 



al = a.i\ . . . a n \ 



C(a,[3) = 



( \ 



yfly 



I a ^ 



where j3 < a. 

We derive a system of equations for 



D a 7t(x), D a n{x). 



We already know that if a is the i unit vector 

D°tt(x) = ir{x) 
D°k(x) = k(x) 

Let the Cauchy data D a n(x) and D a n(x) be known for a = (0, a 2 , a 3 , . . . , a n ) 
and fi(x, k(x)) ^ 0. 

Assume that we have derived algebraic equations D^nix) for D^n^x) and D^k(x) 
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for j3 < a. We apply D a to ( pUPP to obtain 

= ^-(D a 7r(x))f(x,K(x)) (5.1.28) 
+ C{u^){D a -^{x))D?f{x^{x)) 

0<f3<a 

+D a l(x,K(x)). 

This yields an equation for d g™* (x) because fi(x,n(x)) ^ and all the other 
terms d Q° n (x) for i ^ 1 are known from the Cauchy data. 
We apply D a to (|5.1.iq) to obtain 

dx 



= £ C{a^){D a -^{x))D^ gi {x) 



0</3<a 

+D a h(x) + (D a K(x))%(x) 

C(a,f3)(D a ~ l3 K{x))'D l3 l 2 {x) (5.1.29) 

0</3<a 

Notice that this equation ( |5.1.29|) only contains D a n(x) in one term multiplied 
by an invertible matrix so we can express D a function of Z) /3 k(x) for (3 < a 

and -D 7 7r(x). 

In this way we obtain the approximations 

tt(x) « Yl —D a 7r{x)(x - x)) a (5.1.30) 

d \a\<d a ' 

k(x) « ^M^-^r (5.1.31) 

where x is close to x. 



d \a\<d 



5.2 An Example 

Consider the optimal control problem: 

POO 

min / ln 2 (x + 1) + u 2 dt 
u Jo 
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subject to (5.2.32) 
X = xu + u 

x(0) = Xq 



where the domain T> = {x G WL\x > — 1)}. Given the problem ( 5.2.32j ) and the 



schemes described above, we solve the optimal control k(x) and optimal cost tt(x) up 
to degree d and d + 1, respectively; i.e. 

n (x) = n^(x) + n^(x) + ... + ^ d+1 \x) 
k(x) = k^(x) + iS(x) + ... + n [d \x) 

We fix the degree at d = 3 for this example. The interval is T> = (—1,4]. Note that 
the analytic solutions are = — ln(x + 1) and tt*(x) = ln 2 (x + 1). 

5.2.1 Approximation I: Albrecht's Method 



We implement Albrecht's method by using the MATLAB code hjb.m in [T9[ to find 
the coefficients of it and k. We then set-up the polynomials. At each point Xj G V, we 
assign the polynomial approximations 7r 0j = vr (xj) and Ko j = Ko(xj). See Fig.( |5.2| ). 
We solve the optimization problem ( |5.1.23| ). We use the MATLAB code, fmincon.m, 
iteratively to find a point on the level set ir d+1 ^(x) < c; i.e., x G T>i. For this example, 
we march along the x-axis in the both directions. For the interval [0, 4], we march on 
the x-axis towards oo, while on interval (—1, 0] we move towards -1. See Fig.( |5.3|) for 
the psuedocode of the algorithm DRIVER! .m and the actual codes in the appendix. 
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7t(x) of degree d+1=4 




-1 -0.5 0.5 1 1.5 2 2.5 3 3.5 



k(x) of degree d=3 



y 




Figure 5.2: The approximated solutions, n (x) and k (x), via Al'brccht's method compared to the 
real solutions, 7r* and 
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Algorithm DRIVER1 



Input : 



[a,b] , m 



'/interval, mesh size 



f ,1 



'/dynamics, cost functions 



Output: pi(x_j), k(x_j) 



'/optimal cost , control 



Begin 

1. for j=l:mesh 

x(j)=a+ [(b-a)Anesh] *j 
end 

2. initialize f , 1 

4. run hjb.m '/get coefficients 

5. generate pi(x_j), k(x_j) '/make polynomial functions 

6. run endpoint.m, find xbar '/new point on level set 



End 
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5.2.2 Approximation II 

The scheme in the second part is to improve the smooth solutions of Albrecht. See 
Fig.( p.4[ ). The analogue of hjb.m is kovalesky.m in DRIVER2.m, the codes generat- 
ing the coefficients of the polynomials. The coefficients of the polynomials are the 
derivatives and higher order derivatives of 7r and k evaluated at x. First, we solve 
the equations derived from the problem (|5.2.32|) that is described in Section 5.1.3. 
We obtain the equations for the derivatives and its higher order derivatives by using 
MAPLE. Consequently, we write these equations in MATLAB to get the coefficients 
evaluated at x. Then, we piece together the old and the new solutions 7iVi and tci 
and Ki-x and k\ on the interval T>. We denote 7T; as the Ith approximated solution 
and n new as the union of the truncated tti for / = 0, 1, . . . on the interval 
We see that solutions overlapped on some interval around x; i.e. 7T/_i(x) = tti{x) for 
x G [xi — e,Xi] in the case x —>■ oo. In this example, the solutions 7T/_i and ni and 
/t;_i and K\ coincide on some small interval around x. However, we do not expect the 
same result for other examples or in general. In these cases, we take 

ni(x) = min{7ivi(a;),7i-i(x)} 

X 

for some point x on some interval X C T> since the lower approximation is the optimal 
one. The process is then repeated at the next Xi+i. In Fig.( [5.5|) , we compare the new 
approximations, TT new and K new , with the real solutions, 7r* and and Albrecht's 
solutions, 7r and k . The new polynomial, 7r new , is the union of ttq on [—0.6094,0] 
and 71 - ! on [—1,-0.6094]. Similarly, n new is the outcome of glueing k an d «i- See 
Figs. (|irq) and Q 
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Algorithm DRIVER2 

Input : piO , kO '/.Albrecht 3 s approximations 

xbar '/center point of new polynomials 

Output : pi (x_j ) , k(x_ j) '/.optimal cost , control 

xbar '/new polynomials centered around 

Note: We keep the same mesh on [a,b] with m=256 



Begin 

1. for j=l:mesh 

x(j)=a + [(b-a)/mesh] *j 
end 

2. pi~0, k~0, xbar from DRIVERl.m 

3. for 1=1,2,3, . . . 

3a. run kovalesky.m 

3b. generate pi_n(x_j), k_n(x_j) 

3c. if x_j <= xbar 

pi~l(x_j)=pi~{l-l}(x_j) 
k*l(x_j)=pi*{l-l}(x_j) 
else 

pi"l (x_ j ) =pi_n (x_ j ) 
k~l(x_j)=k_n(x_j) 
end °/ if loop 
3d. run endpoint.m, get nxbar 
3e . xbar=nxbar 
end '/for loop 
End 



'/get coefficients 

'/new polynomial functions 

'/glueing solutions if marching — > infinity 



'/newpoint on the level set 
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Figure 5.5: After an iteration on (-1,0], we update n new with 7r_i on (-1,-0.6094]. 
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Figure 5.6: In this case I = 2. Here x\ — 0.5469 and x 2 — 1.3750. 
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7t(x) of degree d+1=4 



20 r 
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x 1 =0.5469 


x 2 =1.3750 


x 3 =2.1875 


x =2.7812 

4 




5 


*0 



























i 


i 


i i 


i 


i 


i i 



0.5 1 1.5 2 2.5 3 3.5 



x 



k(x) of degree d=3 



y 





K 1 


K 2 


K 3 




- - K.M 

K (X) 
new 1 ' 






I 


i i 


i 


i i i 



0.5 1 1.5 2 2.5 3 3.5 4 



Figure 5.7: We iterate up to I = 4 on [0,4]. 
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Figure 5.8: We iterate up to I = 3 on (—1,0]. 



Appendix A 



Codes 



A.l Driver 1 

driver l.m 

7, generates polynomials around via albrecht 
7. 

7, Prager Example 

7, Dynamics: xdot=xu + u; x(0)=xO 

7. Cost: (In x+l)~2 + u~2 

7. Taylorized Cost: x~2 -x~3+(ll/12)x~4- . . . +u~2 

% 

'/, Dimension Guide 

7. d=l,n=l,m=l: input f(l,2), 1(1,3) (linear, starts with quadratic) 

7c output ka=(l,l), py=(l,l) (linear, st w quadratic) 

7. d=2,n=l,m=l: input f(l,5), 1(1,7) (up to 2nd, up to 3rd) 

7. output ka=(l,2), py=(l,2) (up to 2nd, up to 3rd) 

7. d=3,n=l,m=l: input f(l,9), 1(1,12) (up to 3rd, up to 4th) 

7. output ka=(l,3), py=(l,3) (up to 3rd, up to 4th) 

7. d=4,n=l,m=l: input f(l,14), 1(1,18) (up to 4th, up to 5th) 
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output ka=(l,4), py=(l,4) (up to 4th, up to 5th) 



"/Subroutines: 



'/. 



polyl .m 



■/. 



endpoint .m 



'/. 



theG . m 



theF . m 



"/Find coeffients of Taylor''s Expansion 

7,d — degree up to 

y.xbar — Taylor expand around xbar 

function driverl(d) 

mesh=256; 

[xl ,pie ,U,doml ,dom2,rrealu,rrealpy] =polyl (d) ; 
[root] =endpoint (pie) ; 

subplot (2, 1,1) , plot (xl, pie, 'm') ; 
hold; 

plot(xl,rrealpy, 'c — ') ; 

legend( >\pi~0(x) ' , '\pi~-OKx) ' ) ; 

hold; 

title (>\pi(x) of degree d+l=4>); 
xlabeK ' x' ) ; 
ylabel('\pi') ; 
axis([0 4 -2 20]); 



subplot (2, 1,2) ,plot(xl,U) ; 
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hold; 

plot(xl,rrealu, 'c — ') ; 
legend(>\kappa~0(x) > , > \kappa~{*}(x) ') ; 
hold; 

title (' \kappa(x) of degree d=3'); 
xlabeK ' x' ) ; 
ylabeK ' \kappa' ) ; 
axis([0 4 -10 10]); 

Subroutines 

polyl.m 

make polynomial functions 
function [xl , pie , U , doml , dom2 , rrealu , rrealpy] =poly 1 (d) 

n-1; 

if (d==l) 

f=[0 1]; 

1=[1 1]; 
elseif (d==2) 

f=[0 1 1 0] ; 

1=[1 1 -1 0]; 
elseif (d==3) 

f=[0 1 1 zeros(l,4)] ; 

1=[1 1-10 11/12 zeros(l,4) ]; 
else 

f=[0 10 10 zeros(l,9)] ; 

1=[ 1 1 -1 11/12 zeros(l,4) -5/6 zeros(l,5)] 
7. 1=[ 1 1 -1 -1 zeros (1,10) ]; 
end 1 /, ifloop 
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[ka,fk,py,lk]= hjb(f ,1,1,1, d) ; 

"/.generating vectors X, Y (py k ka) 
xslots=nchoosek(n+(l+l) -1,1+1) ; 
yslots=nchoosek(n+(l) -1 , 1) ; 
X=zeros (xslots , 1) ; 
Y=zeros (yslots , 1) ; 

if (d>=2) 
for i=2:d 
xslots=nchoosek(n+(d+l)-l,d+l) ; 
yslots=nchoosek(n+(d)-l,d) ; 
X=[[X] zeros(xslots,l)] ; 
Y=[[Y] zeros (yslots , 1) ] ; 
end'/, dloop 
end /.ifloop 
X=X> ; 

Y =Y> ; 

[xsize , dummy] =size (X) ; 
[ysize , dummy] =size (Y) ; 

aa=0; 

bb=4; 

mesh=256; 

dx= (bb- aa) /me sh ; 

dy=dx ; 

for i=l:mesh 

xl(i)= aa + i*dx; 
a=xl(i) ; 
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X=[a~2; a~3; a~4] ; 

pie(i)=py*X; 

PY=pie(i); 

realpy(i)=xl(i)~2 -xl(i)~3 +(11/12) *xl (i) "4; 

rrealpy(i)=(log(xl(i)+l))~2; 

P=py*X; 

Y=[a; a"2; a~3] ; 
U(i)=ka*Y; 

realu(i)=-(xl(i)-(l/2)*xl(i)~2+(l/3)*xl(i)~3) ; 
rrealu(i)=-log(l+xl(i)) ; 
FK=pragerf(xl(i),U(i)); 
LK=pragercost(xl(i) ,U(i)) ; 
gradP=0 ; 

for k=l:d 

gradP=gradP + (k+1) *py (k) *xl (i) " (k) ; 

gP(i)=gradP; 
end /.gradPloop 

end "/dor loop 

endpoint.m 

'/.finds pt on levelset 
function [root] =endpoint (pie) 

C=fmincon('theG' ,2,0,0,0,0,0,4, 'theF' ) ; 

"/.finding root 
root=100; 
m=100; 
mesh=256; 
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for i=l:mesh 

if (abs(C-pie(i))==0) 

root=i ; 

end 

end '/for loop 

theG.m 

[ka , py] = j ustgivemepk ; 
G=-py* [x~2; x~3; x~4] ; 

theF. m 

y.contraint func for fmincon 
function [F,eqF]=theF(x) 

[ka , py] = j ustgivemepk ; 
u=ka* [x; x~2; x~3] ; 
FK=pragerf (x,u) ; 
LK=pragercost(x,u) ; 
gradP=0 ; 
d=3; 

for k=l:d 

gradP=gradP + (k+1) *py (k) *x~ (k) ; 
end '/.gradPloop 

eps=2~(-6) ; 

Fl=gradP*FK+(l-eps)*LK; 
F2=-(gradP*FK+(l-eps)*LK) ; 
F=[F1;F2] ; 
eqF=0 ; 
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justgivepk.m 

"/output coefficients 
function [ka.py] =justgivemepk 

d=3; 
n=l; 

if (d==l) 

f=[0 1]; 

1=[1 1]; 
elseif (d==2) 

f=[0 1 1 0] ; 

1=[1 1 -1 0] ; 
elseif (d==3) 

f=[0 1 1 zeros(l,4)] ; 

1=[1 1-10 11/12 zeros(l,4) ]; 
else 
7.d=4 

f=[0 10 10 zeros(l,9)] ; 

1=[ 1 1 -1 11/12 zeros(l,4) -5/6 zeros(l , 5)] ; 
end 1 /, ifloop 

7.hjb(f ,l,n,m,d,f_,n_,m_) 
[ka,fk,py,lk]= hjb(f ,1,1,1, d) ; 

A. 2 driver 2 

Driv er2.m 

7. generate polynomial around the x_0 
°/. improves albrecht approx 

y. 
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7, Prager Example 

7. Dynamics: xdot=xu + u; x(0)=xO 

7. Cost: (In x+l)~2 + u~2 

7. Taylorized Cost: x~2 -x~3+(7/12)x~4- . . . +u~2 

% 

7o Dimension Guide 

7o d=l,n=l,m=l: input f(l,2), 1(1,3) (linear, starts with quadratic) 

7. output ka=(l,l), py=(l,l) (linear, st w quadratic) 

7. d=2,n=l,m=l: input f(l,5), 1(1,7) (up to 2nd, up to 3rd) 

7. output ka=(l,2), py=(l,2) (up to 2nd, up to 3rd) 

7. d=3,n=l,m=l: input f(l,9), 1(1,12) (up to 3rd, up to 4th) 

7. output ka=(l,3), py=(l,3) (up to 3rd, up to 4th) 

7. d=4,n=l,m=l: input f(l,14), 1(1,18) (up to 4th, up to 5th) 

7. output ka=(l,4), py=(l,4) (up to 4th, up to 5th) 

t 

7oSubroutines : 

"/, kovalesky.m 

7. poly2.m 

7. glue . m 

7. 

7od — degree up to 

7.xbar — Taylor expand around xbar 

function driver2(d) 

mesh=256 ; 

[xl,pie,U,gP]=polyl(d) ; 



for j=l:4 
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root=endpoint(pie) ; 
ii=root; 

[ka , py] =kovalesky (d , xbar , kl , pil ,pi2) ; 

[xl ,pie2 ,U2,rrealu,rrealpy ,gP] =poly2(d,ka,py ,xbar ,xl) ; 
[newp,newu] =glue(pie2,U2 ,pie ,U, istar ,xl ,rrealu,rrealpy) ; 

figure ; 

subplot (2, 1,1) ; 
hold; 

plot(xl,rrealpy, 'c — ') ; 
plot(xl,newp, 'k') ; 
hold; 

legend ( ' \pi_{*} (x) ' , ' \pi_{new} (x) ' ) ; 
title (>\pi(x) of degree d+l=4>); 
xlabel('x') ; 
ylabeK ' \pi ' ) ; 
axis([0 4 -2 20]); 

subplot (2, 1,2) ; 
hold; 

plot(xl,rrealu, 'c — ') ; 
plot(xl,newu, 'k') ; 
hold; 

legend('\kappa_-OXx) ' , '\kappa_{new>(x) ') ; 
title C\kappa(x) of degree d=3'); 
xlabeK ' x' ) ; 
ylabeK ' \kappa' ) ; 
axis([0 4 -10 10]); 

pie=newp; 
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U=newu ; 

end '/.for 

Subroutines 

kovalesky.m 

"/generates coefficients 

function [ka , py] =kovalesky (d , xbar , kl , pil , pi2) 

y.py — vector containing Taylor coefficients of py 
°/ka — vector containing Taylor coefficients of ka 

ka=zeros(l ,d+l) ; 
py=zeros(l ,d+2) ; 
°/degree=d; 

if (d >= 1) 

[kal ,pyl]=Coeff dl (xbar ,ka,py ,d,kl , pil ,pi2) ; 
end 1 / if loop 

if (d >= 2) 

[ka2 , py2] =Coef f d2 (xbar , kal ,py 1) ; 
end 1 /, if loop 

if (d >= 3) 

[ka.py] =Coeffd3(xbar,ka2,py2) ; 
end 1 /, if loop 

poly 2. m 

"/set up polynomial functions 

function [xl , pie ,U , rrealu , rrealpy , gP] =poly2 (d , ka , py , xbar , xl) 
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mesh=256 ; 

for i=l:mesh 
a=xl(i) ; 

X=[l; (a-xbar) ; (l/factorial(2))*(a-xbar) ~2; (1/f actorial (3) )* (a-xbar) ~3 ; (l/factorial(4))*(a-xbar) ~4] ; 

pie(i)=py*X; 

PY=pie(i); 

rrealpy(i)=(log(xl(i)+l))~2; 

Y=[l; (a-xbar) ; (1/f actorial(2) ) *(a-xbar) '2 ; (1/f actorial(3) )* (a-xbar) ~3] ; 
U(i)=ka*Y; 

rrealu(i)=-log(l+xl(i)) ; 
FK=pragerf(xl(i),U(i)); 
LK=pragercost(xl(i) ,U(i)) ; 
gradP=0 ; 

for k=l:d+l 

gradP=gradP + (k) * (1/f actorial(k) ) *py (k+1) *(xl (i)-xbar) " (k-1) ; 
gP(i)=gradP; 
end '/.gradPloop 

end /.for loop i 

glue.m 

7, attach new polynomial in appropriate interval 

function [newp , newu] =glue (p , u , pie , U, istar , xl , rrealu , rrealpy) 

mesh=256; 

[dum n]=size(p); 

ii=istar ; 
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for i=l:ii-l 

newp(i)=pie(i) ; 

newu(i)=U(i) ; 
end 

newp(ii :n)=p(ii :n) ; 
newu(ii :n)=u(ii :n) ; 

coeffdl.m 

"/.Calculate the deg=l coefficients of Prager''s Example 
function [k,p] =Coef f dl (a,k,p,deg,kl ,pil ,pi2) 

[PP , KK] =Tconst ( a , deg) ; 

p(l)=pil; 

k(l)=kl; 

p(2)=pi2; 

l=(log(a+l))~2 + k(l)~2; 
f=a*k(l) + k(l) ; 

p(3)=(-l/f)*( p(2)*( k(2) + k(l) +a*k(2) ) +2*log(a+l)/(a+l) +2*k(l) *k(2) ) ; 
k(2)= (l/2)*(- p(2) - p(3) - a*p(3)); 

coeffd2.m 

"/.Calculate the deg=2 coefficients of Prager''s Example 
function [k,p]=Coeffd2(a,k,p) 

l=pragercost(a,k(l)) ; 
f=pragerf (a,k(l)) ; 

p(4)=(-l/f)*( 2*p(3)*(k(2)+ k(l) + a*k(2)) + p(2) * (k(3)+2*k(2)+a*k(3) ) +2/(a+l)~2 -2*log(a+l)/ (a+1) ~2 + 2*k 
k(3)=(-l/2)*(p(3) + p(4) +a*p(4)); 

coeffS.m 
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"/Calculate the deg=3 coefficients of Prager''s Example 
function [k,p]=Coeffd3(a,k,p) 

l=pragercost(a,k(l)) ; 
f =pragerf (a,k(l) ) ; 

p(5)=(-l/f)*( 3*p(4)*(k(2) + k(l) + a*k(2)) + 3*p(3)*(k(3) + 2*k(2) + a*k(3)) 

+ p(2)*(k(4)+3*k(3)+a*k(4)) -6/(a+l)~3 + 4*log(a+l)/(a+l) "3 + 6*k(2)*k(3) +2*k(l) *k(4) ) ; 

k(4) = (-l/2)*(3*p(4) + p(5) + a*p(5)); 
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